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CHAPTER ONE 


INTRODUCTION 

The problem of failure In rolling and sliding contacts goes 
back to the period when roller bearings, gear, and railroad 
wheels were first used. Even though many failure modes had been 
analyzed and empirical rules were formulated to guide designers 
against these failures, a new catastrophic failure for bearings 
was recently detected in test rigs for high speed engines. This 
failure mode can be described as a fast growing crack that 
initiates from the surface of an inner raceway of a high speed 
bearing. Fig. 1.1 shows a typical split ring ball bearing used in 
such high speed tests. This type of failure was encountered at 
DN values of three million (where D is the bearing bore in 
millimeters and N is the shaft speed in RPM). These high DN 
values are needed to meet the future requirements of high 
performance aircraft engines. This research investigates this 
failure mode using fracture mechanics theories and numerical 
stress analysis techniques. 

1.1 Failure in Bearings 

Rolling contact devices such as rolling bearings may have 
various causes for failure. The most common causes are plastic 
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indentation, heat imbalance, wear and rolling contact fatigue, 
see Ref . [ 1 ] . 

Plastic Indentation : 

The raceway under a rolling or even under a stationary 
contact undergoes plastic deformation when solid debris is caught 
in the contact region or when the system is overloaded* This 
case produces surface defects that can initiate a fatigue crack 
failure. Fig* 1.2. a shows some typical plastic indentations 
grooves on inner raceway of bearings* 

Heat Imbalance : 

Since all loaded moving contacts generate heat, and on some 
occasion heat is brought into a bearing from the outside, this 
can evaporate the lubricant and soften the contact material, 
which can lead to the gross failure of the bearing. Thermal 
energy must be removed by cooling the bearing to avoid excessive 
temperatures . 

Wear : 

All machined surfaces have asperities that are large 
compared with molecular dimensions* When two surfaces slide over 
each other, the softer asperities either fracture or deform* The 
rate at which these asperities are removed is known as the wear 
rate which usually depend on t lie initial surface roughness, the 
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applied load, the lubricant and the mechanical properties of the 
asperities. Typical surface wear is seen in Fig. 1.2.b. 

Contact Fatigue (Spalling) : 

It is well known that many metals, among them steel, some 
non-metals (ceramics) react to repeated application of the same 
load by gradually building up a fatigue crack. In this failure 
mode a small crack forms in the subsurface, normally associated 
with stress raiser such as void, non-metal li c inclusion, or 
carbide. The crack after repeated loading propagates outward to 
form a spall. But with advances in bearing materials and 
processing techniques, subsurface initiated fatigue is rarely 
encountered. Instead a surface initiated fatigue spalling is 
more likely, where a small crack develops from a debris dent or 
an oxidation pit and propagates inward. These cracks that 
originate at the surface propagate at an acute angle toward the 
bore of the raceway until they meet deformation bands. Here they 
are diverted back toward the surface in the direction the ball 
travels forming a spall. In Fig. 1.2.c, a typical spall, which 
was initiated around a surface furrow is shown. 

While spalling is undesirable, it is a relatively gradual 
fatigue process that can be detected by vibration monitors, chip 
detectors and other oil system monitors. Consequenly the 
affected components can usually be removed before more serious 
secondary damage is incurred. However, at increased rotational 
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speed, (a high DN value), rapid fracture occurs prior to 
significant spalling, and consequently no warning is received 
before the complete failure of the bearing, causing extensive 
damage to the engine* At high DN values, a crack will develop 
from the surface defect of the inner raceway and propagate inward 
at a rate higher than the spalling rate and to a greater depth 
than regular spalling, causing a catastrophic failure. 

1.2 Literature Review 

A large number of investigations have been carried out to 
clarify the mechanism of rolling contact fatigue. However, 
contradiction still persist on the driving force that causes a 
contact fatigue crack. 

From the metallurgical viewpoint it may be considered that 
the process of fatigue consists of two stages; crack initiation 
and crack propagation. The propagation stage has been explained 
so far by two different mechanism. The first mechanism takes 
into consideration the pressure of the lubricant seepage into the 
crack. This mechanism is based on experiments presented by Way 
in 1935, Ref. [2], where he demonstrated that lubricants must be 
present if fatigue crack propagation is to ta*e place. The 
second mechanism is just due to surface traction, where no 
lubricant is present. Fleming and Suh, Ref. 13], were the first 
to analyze crack propagation in sliding contact, under an 
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asperity contact. The tangential hoop stresses due to the 
centrifugal load were completely ignored. Linear fracture 
mechanics was applied using approximate solution for the stress 
intensity factors. Recently, Rosenfield [4] did the same 
analysis with an added feature. He included in his model the 
effect of the friction on the two crack faces under shear, due to 
the high compressive normal stresses. Keer and his coworkers , 
15] [6] [7], analysed the surface crack as well as the subsurface 
crack, under general Hertzian loadings. They also determined the 
crack propagation angles using the maximum tangential stress 
around the crack tip. They determined that the edge crack will 
most likelly grow in a direction of -69 degrees from the 
horizental. In a preceding paper, (6], Keer added extra features 
in modelling the crack propagation of an edge crack inclined at 
65 degrees from the horizontal. In that analysis, the friction 

effect for the cracked faces and the lubricant effect were 

included. A uniform pressure on the open part of the crack faces 
was applied to model the lubricant effect. The variations of the 
stress intensity factors and the angles of maximum tangential 
stresses were described in great detail. But, the investigation 
was not carried further to see what will happen to the angle of 

propagation after the first kink (i.e. crack tip change in 

direction of propagation). All the above investigations were 
based on a two dimensional model of a semi infinite region, and 
not until Murakami et al . , Ref. [8] was a three dimensional model 
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presented- The friction effect on the crack surfaces was 
ignored. It was determined that fatigue crack propagation is 
accelerated predominatly by the hydraulic pressure effect, and 
the direction of crack growth is fixed by the direction of the 
initial crack formed. 

Not until recently were the hoop stresses due to rotational 
speed taken into consideration in the fatigue crack growth. 
Clark, Ref. [9], was the first to suggest that at high rotational 
speed the driving force for crack propagation is the alternating 
mode I stress intensity factor, due to the significant tensile 
hoop stresses. He qualitatively stated that the superposition of 
the Hertzian stress field upon the tensile hoop stress field 
causes an alternating mode I stress intensity factor which leads 
to a rapid crack growth in the raceway. He basically ignored the 
high alternating shear stresses introduced by the Hertzian 
loading as the roller passes from one side of the crack to the 
other . 

Experimental investigations on the endurance and failure 
characteristic of high speed bearings at three million DN values 
were performed by Bamberger, Zaretsky and Signer, [10]. Groups 
of 30, 120 mm bore ball bearings, made of M50 high speed steel, 
were endurance tested at two speeds (12,000 and 25,000 RPM) 
corredponding to 1.44xl0 6 and 3.00xl0 6 DN, respectivally. Initial 
bearing failure at 12,000 RPM was by classical surface rolling 
element fatigue of the inner race. However, continuous running 
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Fig. 1.3 Typical Step-Like Surface Fracture of a Ball Bearing after Failure for a DN Value 
of Three Million, Ref. [51] 
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after a small surface crack had occurred at 3.0X10 DN resulted 
in a destructive fracture of the bearing because of the high hoop 
stresses present due to centrifugal forces. Closer investigation 
of the fracture surfaces showed that the fracture had several 
distinct steps. The initial fracture was normal to the surface 
until it reached the first plateau where it changed direction and 
propagated parallel to the surface in the direction of the ball 
travel. It then, once more, changed direction and propagated 
radially inward until it changed direction for the third time. 
This pattern was observed on several bearings tested, Fig. 1.3. 
This experimental evidence basically contradicts the earlier 
conclusions of Murakami assuming that the crack growth direction 
is fixed by the initial crack direction. 

With the above experimental evidence, one is forced to 
analyze the fatigue crack propagation under the combined tan- 
gential hoop stress due to rotation and the Hertzian load that 
introduces a high alternating shear stresses, as well as high 
compressive stresses as the roller crosses the crack plane. This 
loading condition introduces a complex stress state that leads to 
a mixed mode stress intensity factor. But very little 
information on mixed-mode fatigue crack growth under mixed 
loading conditions is available. The early research on the crack 
propagation direction and time to failure was confined to one 
dimensional loadings. In 1980, Sih, Ref. [11], proposed an 
expression for the mixed mode fatigue according to his strain 
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energy density factor theory: 
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If one considers the crack propagation rate versus the cylic 
stress intensity factor range in a doub]e logarithmic 
presentation, three well known regions of crack growth rate are 
observed, Fig. 1.4, [12]. If the stress intensity range is below 
the threshold range, A K crack growth is not observed. In 
region II, the rate of crack propagation varies lirearly with the 


Lor Crack Growth Rate, da/dN 
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Fig. 1.4 Crack Growth Rate versus Stress Inter.;::;- 
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stress intensity range on a logarithmic scale. In region III, 
the crack growth rate reaches a very high speed since the maximum 
stress intensity factor is approaching the critical stress 
intensity factor for unstable crack growth K^. The most general 
analytical equation for mode I fatigue crack growth rate is due 
to Forman and coworkers [13], where they took into consideration 
all the ranges of the experimental curve, refer to Fig. 1.4: 


da 

dN 
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where R = K T . /K T 

limn lmax 


There are only a few papers in the literature dealing with 
stable crack propagation under mixed loading. One reason is that 
crack growth in the second region apparently occurs under mode I 
conditions, (Fisher [12]), and the main regime for mixed mode 
occurence is region I. This region is the most important for 
fatigue crack initiation and propagation in rolling contact 
fatigue. This may explain why Murakami had to include the 
influence of the lubricant to increase the stress intensity range 
above the threshold range for mode I. 

1.3 Dissertation Outline 

The aim of this study is to analyze the causes that lead to 
the brittle failure of high speed bearings. The object is to 
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provide bearing designers with the tools necessary to cross this 
new barrier for future engines operating at DN values of three 
million and above. 

In order to accomplish this aim, the loadings that are 
significant for the creation of the fast growing cracks should be 
quantified. Aside from the obvious normal Hertzian load across 
the contact interfaces, surface parallel forces also arise from 
friction at the contacts. Also, the centrifugal forces of the 
inner raceway cannot be ignored since they can play an important 
role in the propagation of fast cracks. 

With the external loading being defined, the stresses in the 
inner raceway of a bearing will be calculated numerically using 
the boundary integral equation method (BIEM). The first 
objective of this research is to describe the interaction between 
the Hertz stress field and the centrifugal hoop stress in the 
presence of a surface defect like a pit or a crack in terms of 
stress intensity factors, and The parameters which 
determine the stress fields generated by the rolling contact 
loads are the elastic constants of the bearing materials, the 
design geometries of the raceway and the roller , the 
microgeometry of the surface roughness, the film forming 
parameters of the lubricant, and the geometry of the crack. 

The second objective is to determine the crack growth rate 
and its orientation under this complex loading condition. To 
accomplish this, the major driving forces for fast growing cracks 
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should be identified. It is assumed that combined alternating 
mode I and mode II loadings are behind the step like crack growth 
observed at 3 million DN value. The alternating mode I (crack 
faces moving normal to the crack plane) is associated primarily 
with the oscillation of the high tangential stress due to the 
centrifugal loads from each passage of a heavily compressive 
loaded roller. The alternating mode II (crack faces moving 
parallel to the crack plane) is associated with the change in the 
direction of the shear stress when the roller moves from one side 
of the crack to the other. 

In Chapter Two, the different mixed-mode crack growth 

criteria, for static and dynamic loadings, are presented and 

compared. The crack extension forces are divided into two parts: 

first the region for shear mode growth governed by the shear 

crack extension forces, G . second the tensile mode growth 

r 0 

region governed by the tangential crack extension forces, G 

00 

This crack extension criteria will be used in the analysis of the 
mixed-mode failure of the inner raceway of high speed bearings. 


Chapter Three is 

dedi cated 

to the 

derivation and 

the 

implementation of the 

boundary 

integral 

equation method 

in 


fracture mechanics. 

In chapter Four, the loadings applied on the inner raceway 
of a high speed bearing are quantified. The tangential stresses 
due to the rotational speed and press fit are calculated as a 
function of the DN value. The Hertzian pressure distribution are 
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estimated with and without lubrication. 

Chapter Five gives the variation of the mixed-mode stress 
Intensity factors as function of the Hertzian load, crack length, 
and roller position. Crack extension forces along the crack tip 
are presented as well as the angles of crack propagation, which 
are determined along the direction of the maximum change of the 
crack extension forces. Finally, the life to failure of the 
inner raceway is estimated from crack propagation rate data of 
M50 bearing steel. 

Chapter Six reviews the results and conclusions of the study 
and outlines needs for future research on the mixed-mode failure 
in general and the inner raceway of high speed bearings in 


particular. 


CHAPTER TWO 


MIXED-MODE CRACK GROWTH 

The bulk of fracture mechanics work to date has been 
concerned with single-mode loading. Many practical situations 
are mixed-mode, but pure mode I loading is usually assumed in 
order to obtain easy solutions. This assumption often leads to 
unsafe design. Conversely, neglecting the effects of mixed-modes 
may lead the designer to be over-conservative in an attempt to 
compensate for his lack of knowledge. In this chapter different 
mixed-mode crack growth criteria will be presented and compared 
and the appropriate analysis will be used for the crack 
propagation studies of the inner raceway of the high speed engine 
bearings. 

2.1 Linear Elastic Fracture Mechanics 

Linear elastic fracture mechanics theory is based on an 
analytical procedure that relates the stress and displacement 
fields in the vicinity of a crack tip to the applied load, to the 
structure, to the size, shape and orientation of the crack, and 
to material properties. To determine the stress and displacement 
fields in a cracked body, it is important to define the three 
basic modes of crack-surface displacements, see Fie. 2.1. Mode 


17 


18 





I, opening mode, corresponds to normal separation of the crack 
surfaces under the action of the remote applied tensile stresses. 
Mode II, in-plane sliding mode, corresponds to the crack faces 
moving parallel to the crack plane and normal to the crack 
bottom. Mode III, tearing mode is the out-of-plane shear mode. 
Of the three modes, mode I is technically the most important to 
crack propagation in brittle solids, and has the largest fracture 
tests reported on it for different materials. In this 
investigation, the mixed-mode fracture considered is a result of 
mode I and mode II. This is an in-plane loading condition, and 
the out-of-plane loading, mode III, will be ignored. 

The stress and displacement fields for a semi -i nf i ni t e 
region containing a sharp crack, as given in Ref. [14], have the 
following forms, see Fig. 2.2: 

Mode I: 
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Mode II: 
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Where E Is young's modulus, v is poisson' s ratio, and 

K = (3-V)/(i + V) f or plane stress 

K = (3-4 v ) for plane strain 

Kj and are defined as the stress intensity factors, for mode 

I and mode II loadings respectively. 

The above field equations show that the stress and 
displacement fields in the vicinity of the crack tip are 
functions only of the stress intensity factors and material 
properties. Consequently, the stress intensity factors are 
parameters that represent a measure of the applied stress, crack 
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shape and size, and structural geometry. One should be made 
aware that the above solutions are valid only for snail distances 
from the crack tip. For larger distances from the crack tip 
higher order terms should be included to match the outer boundary 
conditions* 

Since the principle of superposition applies to linear 
stress and displacement fields, Eq . 2.1 and 2.2 are added 

together to determine the mixed-mode stress and displacement 
fields. 


As seen from the above stress fields, the stresses tend to 
infinity as r goes to zero, but the tangential stress and the 
shear stress multiplied by the square root of r go to constant 
values Kj and respectively, f or 0 = 0: 


lim /27rr = K T 

* cu -L 

r •+ 0 

lim a Q /2~Trr = K_ 
. r0 I 

r + 0 


2.3 


As a result the stress intensity factors are measures of the 
intensity of the elastic stress fields ahead of the crack tip. 
Therefore, the stress intensity factors, and , quantify the 
intensity of the stresses for brittle materials, such as the high 
strength steel used in bearings, under mixed-mode loading 


condi tions. 
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2,2 Crack Extension Force 

The crack extension force, also known as the ’’Energy Release 
Rate”, G, is basically derived from the change in the work of 
external loads and the strain energy as the crack grows by an 
amount da, see Fig, 2,3. 


G *= d (- U_ + W T ) / d a 2.4 

t, L 

where is the elastic energy and is the work performed by 
the external forces. This definition is based on the Griffith 
Energy criterion for fracture. 

The above criterion assumes self-similar crack growth, the 
crack extends in the same direction as the main crack* For 
mixed-mode loadings, the possibility for the crack to change 
direction is more likely, as seen in most mixed-mode experiments 
[18], [19], Therefore, the formulation for crack extension under 
mixed-mode loadings, should be able to measure the crack 
extension force at any angle 0 around the crack tip. An 
alternative definition is shown next based on the mechanical 
energy required to close a small crack tip extension along an 
angle 0 measured from the original direction of the main crack. 

The crack extension force, G, can be derived from the energy 
required to open an increment CC* of the crack tip, see Fig. 2.4. 
This definition is taken from Ref. [14] for a self-similar crack, 
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but here it is extended for any kinked crack, to account for the 
possibility that the crack can change direction under mixed-mode 
conditions. The crack extension force is equal to the work that 
would be done by the forces along CC 1 to close the crack 
increment along any angle 6. 

Therefore : 

6 

G = lim 2/6 j h{(a (r) u (6-r)+o (r) u (6-r) • dr 

6 -+ o 0 00 rO r 2. 5 


The 2 is necessary because of the upper and lower surfaces. The 
1/2 accounts for the proportionality between tractions and 
displacements. In this equation only mode I and mode II are 
considered; mode III can be easily added as 1/2 C n u • The 

ZO 2 

relevant stresses to be considered from Eq. 2,l.a and 2. 2. a, are 
those across CC 1 prior to crack extension but in opposite 
direction since a negative force is required to close the crack 
extension CC 1 . 
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The displacements are those across C’C prior to crack closure 


and are based on the new geometry of the crack. 


As an 
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approximation, It is assumed that the displacements for the 
kinked crack have the same general solution as for a straight 
crack, which correspond to the new crack tip C’; therefore the 
displacements for the kinked crack have the following form, based 
on relations in Eq . 2.1.b and 2.2.b for 6 equal tot: 

U r = K II /(2E)( ^ T ^ ( 1+v ) (2 k+2) 

u 0 = Kj / (2E) ( )' 5 (1 +v)(2k+2) 2.1. a 


where and are the new stress intensity factors for the 

kinked crack. 

To find and above, a large number of investigators 

tried to determine the stress intensity factors for a kinked 
crack as function of the stress intensity factors of the main 
crack with no kink. Some of those investigations include the 
works of Russian, Pu and Underwood [15], Wu [16], and Hayashi and 
Neraat-Nasser [17], Hussian et al., Ref. [15], found a 
discountinui ty in the solution as the kinked crack length goes to 
zero. Wu [16] and Hayashi f s [17] results were in agreement 
with each other and they concluded that by taking the limiting 
case, where CC 1 goes to zero, the corresponding relation for 

and K , II are as follows, (see Ref. [17]): 
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2.7.b 


Substituting Eq . 2.6 and 2. 7. a into Eq . 2.5, by making use of 

relations 2.7.b, the crack extension force will have the 
following form: 

G = (1+v) Qc+1) | (Kj cos 3 - 3 K it sin cos 2 | ) 2 + 
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following substitution 
d 4> , gives 
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Substituting the value of the integral into Eq . 2.8, one gets the 
crack extension force for mixed-mode loadings at any arbitrary 
angle 0 measured from the original crack direction: 
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Equation 2.9 is divided here into two parts: first, the 

crack extension force due to the tensile opening displacement 
mode, /ugOggdr, called herein the "tangential crack extension 
force": G qq ; and second, the crack extension force due to shear 
opening displacement mode, / u r ° 0 r dr, called herein the "shear 
crack extension force": G^. This separation of the energy 
release rate, G, is based on the fact that there are 
f ractographi cal ly two types of crack growth, i.e. the shear mode 
growth and the tensile mode growth (Ref. [19]), which are 
controlled by two different driving forces, i.e. G and G , 

re ee 

respectively. 

The separation of the two driving forces for the crack 
growth is given in the following equation for the crack extension 

forces, G Q a and G : 
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When 0 is zero, i.e. looking at the extension forces ahead of the 
crack tip in the direction of the main crack, Eq . 2.10 reduces 
to : 
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Taking a closer look at the tangential crack extension force 

due to the tensile opening mode, Gqq , it can be shown that it is 

similar to the maximum tangential stress theory proposed by 

Erdogan and Sih [18] for determining the angle of crack 

propagation under mixed-mode loadings. 

Otsuka et al . [19] ran experiments on fatigue crack growth 

under mixed-mode conditions and observed two different 

f ratographical modes of failures; the shear mode growth region 

(controlled by G ) and the tensile mode growth region 

r 0 

(controlled by G ), as seen in Fig. 2.5* For given stress 
00 

intensity factors, and K , the dominant mode of failure can 
be determined from the crack extension forces. For very low 
values of the stress intensity factors, i.e. low values of and 
Kji, there is no crack growth. The experimental data points of 
no crack growth are shown as hollow circles. For intermediate 
values of and K^, the crack growth is governed by the shear 
crack extension force, G . The experimental data points of 
shear mode growth are shown as triangles. For high values of 

and Kjj, the crack growth is governed by the tensile crack 

extension force, G . The experimental data points of tensile 

00 

mode growth are shown as solid circles. 

Therefore, below a threshold value of the shear crack 
extension force, G there is no crack growth. Above this 
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threshold value, the shear mode growth is observed. Above a 
critical tensile crack extension force, tensile mode 
growth is observed. 

Another method to account for the crack extension force is 
through the path independent integral, known as the J integral. 
As will be seen in the next section, the J integral is defined 
for self-similar crack extension, i .e . crack extension in its 
original direction . 

2.3 The J integral 

The J integral is often used as an alternative form to the 
energy release rate, G, when the nonlinear material ahead of the 
crack tip cannot be ignored. As first defined by Rice [20], the 
J integral has the form: 

J = J ( W n - t u ) ds 2.12 
s 1 11,1 

where s is a path surrounding the crack tip, W is the strain 
energy density, and t^ is the traction conponent on s, and ds is 
an element of s, see Fig. 2.6. Since Rice, theoretical studies 
have been undertaken to investigate the fundamental mathematical 
characteristics of the J integral, and Sternberg [21] has shown 
that the J integral is a vector of path-independent integrals. 


defined as follows: 
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J k ■ / s ( “ v 'i Vk' <ls 


2.13 


This concept of the J integral may show much promise for mixed- 
mode loadings. Eq . 2.13 defines the extension forces in three 
different directions. Substituting the solution of the stress 
and displacement fields from 2.1 and 2.2 into the J integral, one 
gets, for mode I and mode II only: 


r j .) 

< 

(1+V)(K+1) 
= 4E 

4 

\ K i + K n 1 




l -2 K, K, J 


2.14 


To determine the crack extension at any arbitrary angle, a 
simple trigonometric relation is used to give: 


J 0 9 = COS e + J 2 sin 0 


J 


ee 


(1+V)(K+1) 

4E 


( (Kj + K^) cos0 
- 2 Kj K sin6 


2.15 


Note that for pure mode I or mode II, (K = 0 or K - 0), the 
maximum crack driving force is in the direction of self-similar 
crack growth. This observation contradicts some basic pure mode 
II tests where the angle of crack propagation was observed at 70 
degrees from the original crack plane. 

Bui, Ref. [22], showed that J 2 defined in F.c . 2.13 is not 
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path Independent and that Eq . 2.13 is incorrect. He proposed a 
new approach to calculate the integral by dividing the stress and 
displacement fields into symmetric and antisymmetric parts: 

= 2 . 16 . a 

where: Jj * / s ^ wI n ! “ u \ j ) ds 2.16.b 

r r „H II II , 

and J XI = i s . W n l “ i U i,l } ds 2.16.c 

where u^* and t^* are the symmetric displacement and traction 
fields and u^ and t^^ are the antisymmetric displacement and 
traction fields. Substituting into Eq . 2.16, the closed form 

solutions of the stress and displacement fields around a crack 
tip, one gets: 

(1+V)(K+1) 

AE 

which correspond to the crack extension force for 9= 0, in Eq. 
2.11. As seen from the above definition of the new approach to 
the J integral, only self-similar crack extension is measured, 
and the J integral is an alternative method to determine the 
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stress intensity factors using far away stresses and 

displacements, and not as a crack extension force for any 
arbitrary angle, 

2.4 Mixed-Mode Crack Propagation Criteria 

Many theories had been proposed to predict the angle of 
initial crack growth under mixed-mode loading conditions. The 
most important criteria are the maximum tangential stress [18], 
the minimum strain energy density factor [23], the J integral 
approach [24], the pure mode 1 theory [25], and the crack 
extension forces, as defined earlier. All the methods predict 
almost the same direction of crack propagation for near mode I 
loadings, but when mode II become dominant each method points to 
a different direction. A brief description of each of the above 
methods is given next with a comparison of the first three and 
the crack extension forces over a wide range of mixed-mode 
loadings. 

Maximum tangential stress: 

This criterion is based on the direction of the maximum 
tangential stress around a crack tip [18], which also corresponds 
to the maximum tangential crack extension force, in Eq . 2.11. 

It can be seen that when in Eq. 2.1. a is maximum, o is equal 

vv ^ r0 

to zero. This criterion is based on the observation that the 
crack will grow in the radial direction from the crack tip and 


J/ 


perpendicular to the maximum tangential stress when the shear 
6tress is zero. The direction of crack propagation based on this 
criterion is given by (26]: 


tan(6/2) = 1/4 


V K n 


+ 1/4 




(K I /K JI ) 


+ 8 


2.18 


Miniaum Strain Energy Release Rate: 

Sih [24] has advocated the Idea that the local strain energy 
density factor, S = r dW/dV, should be the governing quantity for 
the fracture process direction, (where r is the distance from the 
crack tip and dW/dV is the strain energy density). The strain 
energy density factor is given by: 


S - a u + 2 a 12 K lKn + a 22 


a = LL ((1 + COS 0 ) (k- cos 0 )) 
11 16 


2,19 


a * — ((sin 0 (2cos 0 - k + 1)) 

12 16 

a = — (( K+ 1) (1 - cos 0 ) + (1 + cos9)(3cos^“ 1)) 

22 i6 


where P is the shear modulus, equals E/(2(l+v)) and K is as 
defined earlier equals (3-4v) for plane strain condition. 

The crack will propagate in the direction of the minimum 
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strain energy density factor given by: dS/d0 » 0. It can be seen 
that the angle of crack propagation will depend on the poisson’s 
ratio* V # in this theory as opposed to the rest of the criteria. 

J Integral Theory: 

This theory is based on the generalized J integral developed 
by Sternberg [21]* and the crack will propagate in the direction 
of the maximum J as defined in Eq • 2.15. The direction of the 
crack propagation based on the J integral theory is given by: 

2 2 

tan 9 = - 2 K i K I] . / < K x + K n > 2.20 

Pure Mode I Theory: 

The pure mode I criterion, from [25], is based on the idea 
that material separation in the process region ahead of the crack 
tip takes place in pure mode I, the opening mode. When this 
criterion is applied* a number of tentative extensions of the 
crack in different directions have to be studied. The direction 
to be chosen has to fulfill the criterion* K = 0. This 
criterion corresponds, to some extent, to the maximum tangential 
stress, where Cqq is maximum and is zero. The maximum 
tangential stress criterion tries to predict the new stress 
intensity factors for the kinked crack from the original crack 
(crack with no kink), while the pure mode I case actually extends 
the crack in the direction for which equals zero. The pure 

mode I theory is more accurate* but needs a large number of 


trials before the direction of * 0 is known. 

Crack Extension Forces Criterion: 

The crack growth direction based on the crack extension 

forces criterion depends on the direction along which the crack 

extension forces are maximum. When the value of the tangential 

crack extension force is less than a threshold value t , the 

60th* 

shear crack extension force is the crack driving force. The 

direction of the crack propagation will be governed by the 

direction of the maximum shear crack extension force, G ^ • For 

r 0 

values of the tangential crack extension force greater than the 
threshold value, the direction along which the maximum tangential 
crack extension force acts, is the direction of propagation, (see 
Fig. 2.5). 

A comparison of the direction of crack propagation for the 

first three criteria and the crack extension forces criterion is 

given in Fig. 2.7 for different ratios of mixed-mode loadings. 

The arc-tangent of the ratio of the stress intensity factors, 

Kf/Kn, is plotted versus the predicted crack propagation angle 

for the above criteria. For near pure mode I loadings, tan 1 
Ki/Kh greater than 60 degrees, good agreement of the three 

methods is observed, but for near pure mode II loadings a huge 

discrepancy is observed. 

For pure mode II, the maximum tagential stress theory as 
well as the maximum tangential crack extension force criterion 
predict an angle of 70.53 degrees, the minimum strain energy 


Predicted Crack Growth Direction, (deg.) 
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Fig. 2.7 Prediction of the Crack Growth Direction for Different 
Criterions 
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density factor theory gives an angle of 83.61, while the J 
integral approach gives a zero angle. The difference between the 
J integral and the other criteria can he explained from the wrong 
definition as seen earlier in section 2.3. 

For pure mode II, based on the crack extension force 
criteria, the propagation direction will be zero degree if the 
shear crack extension force Is governing and 70.53 degrees if the 
tangential crack extension force is dominant. 

The crack extension force criteria will be used to determine 
the crack growth rate and direction in the failure analysis of 
the inner raceway of high speed bearings. 

Knowing the crack extension forces and the fracture 
criteria, a crack propagation law can be formulated to predict 
the crack propagation rate under mixed-mode loadings. 

2.5 Fatigue Crack Propagation Law for Mixed-Mode 

Under cyclic loading the crack propagation occurs at lower 
stresses than it does for static loading. For typical fatigue 
fracture of brittle materials, macroscopic plastic deformation is 
not involved. Instead, submicroscopic slip lines developed ahead 
of the crack tip, that can intensify and broaden with each cycle. 
Those slip lines are developed along unfavorable oriented grains. 
But these slip bands cannot cause fracture in themselves, since 


their effect is to harden the material • When a material work- 
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hardens, it will be vulnerable to the tangential stresses around 
the crack tip, which try to open up the defects causing the 
increase in microcrackings with increasing number of cycles. The 
formation and the growth of these microscopic cracks ahead of the 
crack tip represent the damage accumulation. This damage 
accumulation is the cause of the weakening of the material ahead 
of the crack tip, and the decrease of the crack resistance force. 
The crack resistance force should be a decreasing function with 
increasing number of cycles, and how fast the crack resistance 
force decays is a function of the maximum range of G or K. But 
the crack growth is not continuous since experimental evidence 
has shown that the crack growth stops or slows down to a level 
where no crack growth is observed. After N numbers of cycles 
when the crack resistance ahead of the crack tip is decreased, 
the crack will grow instantaneously a distance and then stop. 
At this position the crack tip has reached a region where the 
damage accumulation is non-existent or is so low that more cycles 
are required before the crack can grow. This explaine the 
experimental observed jump like growth of fatigue crocks. 

Generally, one can find the formulation of the fatigue crack 
growth as follows, Ref, [27]: 

d a(t) = F ( a(t), P(t), d(t), t) d t 2.21 


where a(t) is the length of the dominating crack, p(t) is the 



stress load and d(t) includes the damage accumulation and 


material properties • 

A simplified version of Eq . 2.21 was given by Paris [28], 
where the damage accumulation was ignored and the load and crack 
length were combined together as a function of the stress 
intensity factor. Thus, Eq . 2.21 reduces to, (28]: 

da / dN = C ( AK )" 2.22 

Eq. 2.22 is based on experimental results, for mode I loading 
condition, and basically describe region II, in Fig. 1.4. But 
more careful experiments had shown that the mean stress has an 
effect on the fatigue crack growth, and Eq . 2.22 was therefore 
modified to the following form, Ref. [29]: 

da / dN = C ( AK ) ni ( K )° 2 2.23 

max 

Recently, Forman et al., Ref. [13], have provided a general 
expression which takes into account all the regions of Fig. 1.4; 

da / dN = C ( AK - A K^ ) ni / ( (1-R) K^ -aK)" 2 2.24 

where R is the stress ratio, and equals to K /K . 

min max 

All the above equations are based on pure mode I loadings, 
and must be generalized to include the mixed-mode loading 
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conditions and the damage accumulation. 

The first mixed-mode fatigue crack growth equation was 
derived by Sih [11], based on his strain energy density factor, 
where fatigue crack growth occurs in the direction of the maximum 
change in the strain energy factor, and the fatigue growth 
equation becomes: 

da / dN = C ( A S , ) n 2.25 

rain 

where A S , is given by [11]: 
min 


AS . = 2(a. . K AK + a. A K AK + K A K ) 

rain 111 I 12 II I I II 


* a 22 K II iK II> 


2.26 


where K is the mean stress 
intensity factor range, the 
crack propagation for which 
Eq. 2.19. 


intensity factor, A is the stress 


a^'s are functions of the angle of 

AS , is a maximum and are given in 
rain 


Alternatively, one can easily extend the crack propagation 
to the mixed-mode loadings, by making use of the crack extension 
forces along any angle 0, instead of the mode I stress intensity 
factor in the equations on page 43. By substituting relations 
2.10, for the crack extension forces in Forman's equation, the 
crack growth rate will be equal to: 

For Tensile Mode of Crack Growth: 
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Therefore, depending on the mixed-mode loading conditions, one 
must choose the proper mode of failure. Eq . 2.27 would he valid 
for the tensile mode of fatigue crack growth while Eq . 2.28 would 
be valid for the shear mode of fatigue crack propagation* The 
is the threshold crack extension below which no crack 
propagation is observed. 

To determine the dominant crack growth mode, the values of 

the maximum changes of the crack extension forces are compared 

with the threshold values. Fig. 2.8 represents the boundary for 

different crack growth modes. For values of AG C _ greater than 

A^QQthtens ion mode is observed and Eq . 2.27 is valid. For values 

of AG less, two regions are observed; a) no crack growth region 

for AG „ less than AG b) shear mode growth for value of AG 

re r 0th e a r6 

greater, and Eq. 2.28 is valid. 

This separation of the growth modes can easily explain the 
different crack growth directions observed in pure mode II tests 
ran by Buzzard et al., Ref. [52]. For very small cracks the 




47 


crack extension forces are so small that the shear mode is 

dominant causing the crack to extend in the original direction. 

For larger cracks, where the tangential crack extension are 

higher than the threshold AG _ f the crack changes direction, 

00 th’ 

and propagates along 70.53 degrees, in the direction of the 
maximum change in Ggg, indicating a tensile growth mode. 

Thus, the crack growth direction will be determined from the 
angle 6 along which the change in the crack extension forces are 
maximum. The crack will keep on growing in this direction until 
a new loading combination changes the direction of the maximum 
change of G^. 

The crack extension forces criteria will be used in 
analyzing the mixed-mode failure of the inner raceway of the high 
speed bearings* But before determining the crack growth 
direction and rate under mixed-mode loadings, the stress and 
displacement fields have to be determined first using the 
boundary integral equation method for any arbitrary crack 
geometry and orientation, and loading conditions as will he seen 
in the next chapter. 


CHAPTER THREE 


Boundary Integral Equation Method 
3.1 Introduction 

Methods of analysis in elasticity, and in most other 
scientific and engineering fields, have been revolutionized with 
the advance of computers. Most of the solutions, two dimensional 
or three dimensional, had been obtained for infinite or semi- 
infinite bodies using the stress functions techniques, which 
satisfy the desired boundary conditions near the origin and have 
the properties that the stress and/or displacement vanish or 
remain bounded as the boundary at infinity is approached. With 
the advance of computers, solutions for finite geometries and 
mixed boundary conditions were attainable by numerical 
techniques. First, the finite difference method was used, by 
replacing the differential equations by their classical finite 
difference equivalent. More recently, the finite element method, 
which makes use of the variational statement of the original 
differential equations to obtain solutions, has gained 

popularity. In both methods, finite difference and finite 
element, the continuum is discretized, making the accuracy of the 
solution a function of the fineness of the discretization. 
Another method of analysis, recently rediscovered by Rizzo in 


48 



49 


1967, Ref. [31], the boundary Integral equations method, BIEM, 
offers an attractive alternative. The boundary integral method 
involves the transformation of the partial differential equations 
describing the behavior of the unknowns inside and on the 
boundary of the domain to integral equations over the boundary, 
i.e. the integrals are functions of the boundary data only; thus 
enabling the reduction of the dimensionality of the problem. 
Closed form solutions can be obtained for simple finite 
geometries and loading conditions, when the integral equations 
are solvable. The resulting solution is the exact solution of 
the differential equation for the given boundary conditions. For 
complicated geometries and loadings conditions, the integral 
cannot be solved analytically and approximations have to be 
introduced. Therefore, in the BIEM, inaccuracies arise from 
numerical integration procedures, which means that by refining 
these approximations any degree of precision is theoretically 
achievable. By the approximation of the integration, the 
integral equation will be transformed to a set of linear 
algebraic equations. The resulting system of equations is 
smaller by an order of magnitude than those for the finite 
difference or the finite element methods, but are fully- 
populated, whereas in other methods the matrices are symmetric 
and most of the time banded. 

The boundary integral equation method has seen increasing 
popularity in recent years because of the some of the advantages 
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listed below, Ref- [32], [33]: 

i) a smaller set of algebraic equations to solve, ii) simple data 
preparation to run the problem, only boundary discretization is 
needed, iii) infinite and semi -i nf ini te problems are properly 
modelled, iv) accurate solution of selective internal stresses 
and displacements, and v) good resolution for stress concen- 
tration problems. 

In view of the above advantages, the boundary integral 
equation method was used herein to analyze the stress and 
displacement fields around a surface crack of the inner raceway 
of the high speed bearing. 

A review of the two dimensional plane strain elastostatic 
solution by the boundary integral equation method is given next, 
with its implimentation in a general computer program. The three 
dimensional derivation is similar to the two dimensional and its 
implimentation can be found in Ref. [34], 

3.2 Mathematical Derivation 

The most direct derivation of the boundary integral 
equations is based on a singular solution of the Navier 
equations. The Navier equations of equilibrium (in terms of 
displacements) for two dimensional problems in elasticity are, 
Ref. [32]: 




52 


V 2 u. + l/(l-2v) 6 - 0 

i > * 

e - UjJ *.J - >.2 


3.1. a 


where are the displacements, v Is Poisson’s ratio and the 
usual tensor notation is used, where a repeated subscript 
indicates summation over its range and a comma indicates partial 
differentiation. The Navier equations can be written in another 
form as : 


u i,jj + 1/(1 " 2v) “k.ki * 0 


3.1.b 


A solution to these equations can be obtained by making use 
of Kelvin's singular solution due to a single unit concentrated 
force acting in the interior of an infinite body [35] (see Fig. 
3.1). 

The displacement field at any point Q at distance r from the 
point P, where the force is applied is given by Ref. [35] : 


Uj - (-1/[8TT6(1-V)] }|(3-4V) lnr 6^ - r ± r^j 3.2 


. a 


or 


u i • e j 


3. 2. b 


where e. are the unit base vectors, and 
J 


l). . * 1-1 /(8 ttG(1-v)] } | ( 3-4v) lnr 6 - r r 

i j 1 1 J » 1 > J 


3.2.c 
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where r Is the distance between a field point Q with coordinates 

Q p 

(Xj ) and the point of load application P with coordinates (x^ ); 


r - ( (x^ - x^) (x^ - Xj P ) 


*5 


3.3 


For plane stress, poisson's ratio v Is replaced by v /(1 + v), in 
all the equations, Ref. [32]. 

If we consider the field point Q to he on the boundary of a 
body cut out of the infinite region, then the traction forces can 
be determined on this boundary by 


t . = a. . n . 
i ij 3 


3.4 


where n^ is the component of the outward norial in the j 
direction at the surface of the body. Expressing the stress 
tensor in terms of displacements [35] 

°1J ‘ { 2Cv/(l-2v)} S tj ., ,*t <“ jfi * "j.j) ’-5 

where G is the shear modulus, differentiating Eq . 3.2 and 

substituting in Eq . 3.5, Eq. 3.4 becomes: 


tj - | — ( 1 — 2 v) / [ A tt ( 1 —v ) r ] J jdr/dn [6 4j +[2/(l-2.)] r ^ r^] 


- r ,i n j + r ,j "i t 6 j 


3. 6. a 


or 


*i • T ij 


3.6.b 
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where 


Tjj - |-(l-2v)/ (4tt ( 1-v )r] J jdr/dn +(2/(l-2v)] r j r jj] 


- r .i "j + r .j "•) 


3.6.c 


We now make use of Betti’s reciprocal theorem [36] which 

states: If an elastic body is subjected to two systems of surface 

* 

tractions and t y then the work that would be done by the 

* 

first system t^ in acting through the displacement u ^ of the 

second system is equal to the work that would be done by the 

* 

second system t j acting through the displacement of the 

first system, i .e. ; 


l 


'j u 3 


J. 




3.7 


where s is the boundary surface of the body, and ds is an element 
of surface area. 

Suppose we now choose the second system of traction and 

* * 

displacement ( t and u ^ ) to be the one produced by a single 

unit concentrated load, and the system u^ , t^ to correspond to 

the solution we are seeking. Since we know the solution to the 

unit concentrated load (Kelvin solution), we can solve for any of 

the unknown traction and displacement (t., u.) by substituting 

* * 

Eq. 3.6 and 3.2 for t . and u , respectively ar.d solving the 

J 3 


integral. Eq. 3.7. Because of the singular nature of U and Tj. 

at r ** 0, it is necessary to employ a limiting process as shown 
in Appendix A, resulting in the following equation, known as the 
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boundary Integral equation: 


f t U e ds = f u. T e ds - C u 
J s jjll •'sjjll Jli 


3. 8. a 


or In another form 


C u = / U t, ds - / T u. ds 

ij j s Ij j s Ij j 


3.8.b 


where C, . = 6 for Interior points and C.. = 6 /2 for boundary 

ij ij ij ij 

points with smooth tangents. Eq • 3.8 is also known as 

Soraigl i ana 1 s identity. For very simple geometries and boundary 


LIUUU-l L i U1I3 JUiUigiiaua o iuc 


II L 1 L ) 


OU L 1 J 1 CIV. L U I 


. J _ J . 

“*•> “£> 


analytical solutions, but for complex bodies a numerical solution 
is necessary and is discussed in the next section. 

Once the unknown traction and displacement are determined on 
the boundary, internal displacements and stresses can be 

calculated as functions of the boundary displacements and 
tractions. For internal displacement, Eq . 3.8 is used with C 

® j , however for internal stresses Eq . 3.8 is differentiated 

and substituted in Eq. 3.5, to get [32]: 


ij 

where 


fs U ijk C k ■ 


T ijk u k ds 


3.9 


U ljk - |(l-2V)/l4Ml-v)r]| |« lk r ^ + ^ r >k 


,i T .j r ,k 


+ 2 r 


/ ( 1-2V) 
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T ijk " |G/[2TT(i-v)r 2 ]) 1 2 dr/dn l(l-2v) + ; Jk r 

+ V r .i - 4 r ,i r .j r ,k' 

+ "i < 2vr ,j r ,k * ( ’- 2v) V 

+ "j ‘ 2vr ,i r ,k + °' 2V) V 

+ "k I2( 1-2V ) r , i j - ( l-4v ) 

Thus, the displacements and stresses at any interior point, if 
needed can be obtained by integrating numerically the boundary 
equations, Eq • 3.8 and Eq . 3.9, respectively with = ^ij* 

from the solutions of the stresses and displacements at the 
boundary . 

3.3 Reduction of the Integral Equation to a set of linear 
Simultaneous Equations 


The first step in solving the boundary integral equations is 
to reduce them to a set of linear simultaneous algebraic 
equations, if the integral are unsolvable in closed form. The 
boundary of a body to be analyzed is divided into M surface 
elements. Those elements can be linear or curved to map the 
boundary geometries, as seen in Fig. 3.2. Eq. 3.8 can then be 
rewritten as: 

M r M 

c ij U j is. u ij 'j ds " - £ L T n “j 3 -‘° 

b=l b=l J ‘ 


As an approximation the traction, 


and displacement , 


u . are 
J 
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b) Linear Variati:r. 


c) Quadratic Variation 


Fig. 3.3 Different Element Types 
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assumed constant or to vary linearly or quadrat icaly over each 
element, as seen in Fig. 3.3. Thus, the tractions and 
displacements in each element can be approximated by: 
tjCx) = N k (x) t 1 

Uj(x) = N^(x) 3.11 

For constant elements: 

k = 1 and N=1 3.12 

Here the traction and displacement are approximated as a constant 
having the values at the mid-point. 

For linear elements, Ref. [37]: 
k - 1,2 

N 1 = (l-x)/2 N 2 = ( l+x)/2 3.13 

J(x) = (x2-xl)/2 

where J(x) is the Jacobian that transfers the integral from the 
global to the local coordinate system, where xl and x2 are the 
in-plane coordinates of the nodes. 

For quadratic elements, Ref. [37]: 
k = 1,2,3 

N 1 = x(x-l)/2 N 2 = x(x+l)/2 N 3 = (x+l)(l-x) 

J(x) - (x2-xl)/2 + x(xl+x2 - 2 x3) 3.14 

If the surface is represented by M elements, then the integral 


equat i ons become : 
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'ij »j(P) +EL u (Q bk ) /t (P,Q bk ) N k (x) J b ( 
b= 1 k=l J ^ 

M n b r 

L L t (Q bk ) / U 
b=l k=l 


x) dx 


(P,Q bk ) N k (x) J b (x) dx 3.15 


where n b Is the total number of nodes In element b. J b is the 
Jacobian function for the bth element. The terms, Uj(Q bk ) or 
tj(Q bk ), are the nodal values of the displacements and tractions, 
respectively, for the kth node within the bth element. Eq. 3.15 
should be repeated N times, corresponding to the total number of 
nodal points. The total integral equations can be written as: 


N r 

v V L , _ a x , _ a 
N M n 

EEL v< bk > 

a=l b=l k=l 


« ' \ ■ ..me. i _ ,_a _nx. 

' _r Z-rZ^ U i W ' I 1 ii vr ' 

b=l k=l 

J ' Ujj(P a ,Q bk ) N k (x) J b (x) dx 


k. . _b. . 

« v x ; j v, x ) 



3.16 


The expression in Eq. 3.16 represents a set of 2N equations which 
can be written in matrix form as: 


+XI abk N k J b dx J |uj a |-Jx z/u ij abk N k J b dxj jt^j 3.17 

where i,j = 1,2, a = 1,N, b = 1,M, and k - l,n b 
or in general form 

[ A ) ju) - r B ] jt) 


3.18 
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The integrals in Eq . 3.17 are easily solved numerically or 

in closed form over the domain of each element. For Q f P , an 

eight point gaussian quadrature formula is irapl i mented in the 

bk a 

computer program, while for Q = P , closed form solution for 
the lnr singularity is calculated. 

For the case of a traction problem where the t’s are known, 
or the case of a displacement problem where the u's are known, 
Eq. 3.18 reduces to the form 

[ A ] (x) = jc) 3. 19 

Eq. 3.19 represents a set of 2N linear algebraic equations which 
are to be solved by Gauss Elimination method. In case of mixed 
boundary value problem, where some values of both t and u are 
specified, it is necessary to interchange the columns of matrices 
A and B (in Eq. 3.18), so that all the unknown quantities are 
contained in the column vector u and the known values are 
contained in t, before reducing the equation to the form of Eq . 
3.19. 


By placing the nodes at the corners of elements two 
difficulties become apparent: 

1) The possibility exists for nodes to be placed at sharp 
edges of the body rather than at flat surfaces. C^, * n ^9 • 
3.16, is equal to 1/2 6 for flat surfaces. For nodes at edge 

discountinuities can be computed in two ways. One uses a 

limiting process as derived in Appendix A. The second method is 
based on rigid body motion, as explained in Ref. [38], and is 
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Fig. 3.4 


Traction Discontinuity Models 
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given below: 


N 

Cj^P.Q) - 

b=l 



(P,Q) 


dx 


for P h 0 b 3.20 


In the computer program, all . terms are computed using Eq . 
3.20. The value for on flat surfaces was computed by the 

above equation and it was found to be exactly equal to 1/2 6^ as 


predicted by the analytical formulation of Appendix A. The 
second method was implimented in the computer program merely as a 
check to the numerical integrations. 

2) Placing nodes at corners of elements assures the 
continuity of displacements and tractions. However, in modelling 
real problems a step change in traction may exist. To assure 
discontinuity of applied tractions, the input values of traction 
are associated with the element they act on instead of the nodes. 
As an example, consider two adjacent elements which lie in two 
different planes ( see Fig. 3.4). Element 1 is under uniform 
tension t while Element 2 is traction free. If the traction is 
associated with node A directly, an extra shearing traction 
exists in Element 2 varying from zero at node C to t at node A. 
By assigning the traction to a node of a specific element, in 
this example to node A of Element 1, the problem of adding extra 
traction is avoided. Alternatively, one can place two distinct 
nodes between elements 1 and 2, Fig. 3.4.b, but this method is 
not implemented here. 
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Some simple examples are shown next to test the boundary 

integral equation computer program written for the CRAY-XMP at 
NASA Lewis Research Center. 

3.4 Uniform Distributed Load 

The computer program of the two dimensional Boundary 

integral equation method (BIEM2), is first applied to a 

rectangular body with uniformly distributed vertical pressure on 
part of the boundary, under plane strain conditions, see Pig* 
3.5. To simulate the serai -i nf i ni te region a uniform reaction 
force was applied to the bottom part of the plate. The 

dimensions of the body are W/b - H/b = 30, where H is the height, 
W is the width and b is half the length of the pressure applied. 
The poisson’s ratio is equal to 0.3, and the applied pressure 
load is normalized with the modulus of elasticity, p/E ~ 1* 

Two different sets of runs are studied. In the first set 
linear elements are used, while in the second set the quadratic 
elements are impl iraented . Different number of nodes are used in 
each set to analyze the order of convergence of the specific 
element type. A mesh generating program was written to 
facilitate the input data for the boundary discri tazation. The 
required inputs are the material properties, the number of major 
subdivisions with the number of nodes desired in each one of 
them, and finally the coordinates of these major subdivisions. 
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The computer program (GENE 1 ) will then generates the coordinates 
of the nodes and the element formations, with the proper fixities 
and loading conditions. Any loading condition can be implimented 
by just changing the loading function in the program. 
Coordinates of internal points, if desired, are also generated. 
For the above example, nine major subdivisions are chosen, (see 
Fig. 3.6), with fixities in the y-direction at coordinates 2 and 
6, and in the x-direction at coordinate A. A uniform traction is 
applied between coordinates 8 and 9, and along the bottom plate 
from coordinate 3 to 5. The numbers of nodes used in each 
subdivision is shown in Fig. 3.6. For the convergence analysis N 
in each mesh was 2, A, 8, 16, 2A, and 32 respectively. 

To study the convergence of each element type, the values of 

the stresses at five internal points, picked below the edge of 

the uniform load, at distances y = b, and x = b, 2b, 3b, Ab, 5b, 

respectively, were compared, (see Fig. 3.6). Table 3.1 shows the 

values of the stresses for different mesh size at x = b and y - 

2b. Good agreement with the closed form solution is observed for 

N = 32, for a and o , but for o the error is much higher 
* xy y x 

because of truncation error. 

A plot of the average error for the stresses ( 0^, a xy ), 

as function of the number of nodes, for the two type of elements 
is shown in Fig. 3.7. The average error of all the stresses for 
those five points is obtained from the elastic solution of a 
semi-infini te region and the analytical results from the B1EM 
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Table 3.1 Variations of the Stresses o , o , and o , as a 

x* v* xy 

Function of the Mesh Sizes, for an internal point 
at x = b and y = 2b of the uniform loaded plate. 


a x /p °y /p C xy /P 


N 

linear 

Quad . 

Li near 

Quad . 

Li near 

Quad . 

2 

-0.1142 

-0.1195 

-0.3718 

-0.3840 

0.1205 

0.1326 

4 

-0. 1088 

-0. 1063 

-0.4020 

-0.4117 

0. 1336 

0. 1412 

8 

-0. 1002 

-0.0980 

-0.4138 

-0.4165 

0.1479 

0.1531 

16 

-0.0972 

-0.0963 

-0.4128 

-0.4124 

0. 1561 

0.1588 

24 

-0.0969 

-0.0964 

-0.4120 

-0.4112 

0. 1573 

0.1593 

32 

-0.0967 

-0.0966 

-0.4105 

-0.4098 

0. 1586 

0.1595 


exact -0.0908 -0.4092 0.1592 
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program. As seen in the figure, the difference between the 
linear and the quadratic elements is very small for the same 
number of nodes, and both methods converge to the same answer. 
This shows that for this problem, the converged solution does not 
depend on the element type. The advantage of the quadratic 
elements over the linear elements comes into consideration in the 
CPU time required to set up the integral equations. Since the 
integrals are evaluated over each element, and the number of 
elements for the quadratic type is half the number of elements 
for the linear one, the time needed is cut in half. 

Fig. 3.8 shows the CPU time in seconds required on the CRAY 
XMP, for the linear and quadratic elements as a function of the 
number of nodes. As seen in the figure the average saving in CPU 
time between the quadratic and linear elements is about 33%. The 
reduction is only 33% because the CPU time reported is the total 
CPU time needed to set up and solve the system of algebraic 
equations, and both element types need the same solving time, for 
the same number of nodes. Considering the saving in CPU time, 
quadratic elements were used in the application of the BIEM to 
the mixed-mode failure of the inner raceway of the high speed 
bearings, and in the following examples on the stress intensity 


factor calculations 
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3,5 Determination of the Stress Intensity Factors 

The next examples relate to the determination of the stress 
intensity factors for bodies with cracks using the BIEM. There 
is extensive literature on the numerical determination of the 
stress intensity factors. One of the methods which was adopted 
from finite elements is the use of a singularity element at the 
crack tip, Ref. [39]. The simplest crack singularity element is 
the quarter point element, where the mid-point of the quadratic 
element at the crack tip is moved to the quater point position, 
Ref# [39]* The disadvantage of these elements is that a 
transition element is needed and the crack tip should be fixed to 
give accurate results for the stress intensity factors. A second 
disadvantage is that the length of the crack tip element affects 
the result. Therefore, the length of the quarter point element 
has to be adjusted until the desired accuracy is achieved. Due 
to these disadvantages, this method is not implemented here. 

The second method to determine the stress intensity factors 
is based on their basic definitions; as the limit of the stress 
multiplied by the square root of the distance to the crack tip, 
for the stress ahead of the crack tip, at a distance of 0.2 of 
the crack length to 0.5 of the crack length, (See Eq. 2.3). In 
this method the stresses ahead of the crack tip, the tangential 
stresses and the shear stresses, are multiplied by the square 
root of 27Tr and plotted versus r. Curves are then fitted 
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/3 


throught those points and the values at r = 0, equal the stress 

intensity factors, and respectively. This method is also 

known as the extrapolation method for determining the stress 

intensi ty factors . 

Finally the J integral as defined in section 2.3 can be 

evaluated along any path around the crack tip, and the values of 
the resulting integrals are proportional to the stress intensity 
factors. The extrapolation method and the J integral method are 
compared next for a single edge crack under mixed-mode loading 
condi t i on . 

The single edge crack under uniform tensile load is analyzed 
first, with dimensions a/H = 1/15 and W/H = 1, where a is the 
crack length see Fig. 3.9. The tensile stress, t, is normalized 
with the modulus of elasticity E, t/E = 1. The number of nodes 
in each subdivision are shown in Fig. 3.9. The number of nodes, 
N, along the crack, is varied between 20 and 140. The mode I 

stress intensity factor, , was calculated using the 

extrapolation method using internal points ahead of the crack tip 
between 0.2a and 0.5a. For points closer than 0.2a, the values 
of the stresses are not accurate. The variation of K^. with the 
number of nodes along the crack front is shown in Fig. 3.10. An 
asymptotic value for the stress intensity factor is reached for 
values of N greater than 100. Comparing the asymptotic value 
with the published solutions, Ref. [40], the difference between 
the two values is less than 1%. 


The J integral was also 
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Fig. 3.11 The Three Paths Considered for 
Calculat ions . 


A 


* 



the J Integral 
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Table 3.2 Variation of the J Integral and the Corresponding 
stress Intensity factor for different paths, for 
the Pure Mode 1 Loading Condition. 


Path 

JjE 

K I 


t 2 a 

t /a 


Eq. 2.13 

Eq. 2 

1 

3.8052 

2.045 

2 

3.8113 

2.046 

3 

3.8146 

2.047 


extrapolation 2.043 

Ref. [40] 


2.032 


calculated for N equals 100. Three different rectangular paths 

were chosen as shown In Fig. 3.11. To determine the values of 

u . in the J integral formulation, neighboring points were also 
i * * 

used to determine numerically the values of the derivative of the 
displacement. The values of J integral for the corresponding 
paths are shown in Table 3.2, based on Eq . 2.13. Also shown is 
the corresponding mode I stress intensity factor using Eq . 2.14, 
and the results using the extrapolation method as well as the 
published solutions from Ref. [40]. As seen from Table 3.2, the 
results of the three different paths agree well with each other 
and also give the same result as the extrapolation method. Thus, 
the J integral and the extrapolation method give the same value 
for the stress intensity factor for this pure mode 1 case. 

The second example is for a plate with a crack under shear 
mode loading, pure mode II, see Fig. 3.12. The dimensions of the 
plate are the same as for the mode I case; a/H = 1/13, W/H = 1. 
Two antisymmetric uniform applied pressure loads are located at 
equal distance from the crack plane. The distance y^ of the 
raid-positions of the pressure load to the crack plane equals 
0.25a. The magnitude of the pressure P applied equals 50E, and 
acts on a segment of length equals 0.01a. The variations of the 
mode II stress intensity factor as a function of the number of 
nodes along the crack front, using the extrapolation method, is 
shown in Fig. 3.13. An asymptotic value for the mode II stress 
intensity factor is reached at values of N greater than 100. The 
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Fig. 3.13 Variation of the Mode IT Stress Intensity Factor as a Function 
of the Number of Nodes N, Along the Crack Front, (p/t=50) 
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J integral method was also calculated using both formulations as 
describe in section 2.3, Rice formulation Eq . 2.13 and Bui 
formulation Eq . 2.16, for N = 100. For Bui’s formulation the 
stresses and the displacements along any path had to be devided 
into symmetric and antisymmetric parts: 


Di spl acement s : 


Stresses : 


r i i 
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ii 
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xy J 


a + 
xy 



Sycmet r i c 


Anti symmet ri c 


Sy-metric 


Anti symmetric 


where the 1 indicates the component of the stress or the 
displacement at the opposite position, with respect to the crack 
plane, of the point under considrerat ion, see Fig. 3.11. Table 
3.3 represents the values of the J integrals for the two 
formulations: Rice (Jj, J 9 ) and Bui (J^, Jjj), along the three 

paths given in Fig. 3.11 for the pure mode II example. 
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Table 3.3 Values of the J Integral for the Three Paths, for 
the Pure Mode II loading Condition, (p/t = 30). 



Ri 

ce 

Bui 

Path 

J j E 

j 2 e 

V : J n E 


t 2 a 

2 

t a 

t 2 a t 2 a 


Eq . 

2. 13 

Eq. 2.16 

1 

0.40877 

0.00 

0.00525 0.40353 

2 

0.44304 

0.00 

0.00168 0.44217 


3 


0.39522 


0.00 


0.00074 


0.39449 
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Table 3*4 Values of The Corresponding Mode II Stress 

Intensity from the J Integral Calculations, and 
the Value from Extrapolation method, (p/t = 30). 


Ri ce Bui 


Path 

K I 

K II 

K I 

K II 


t y/5" 

t Ja 

t v /cT 

t/a 


Eq . 

3.23 

Eq . 3. 

24 

1 

0.00 

0.6702 

0.07596 

0.6659 

2 

0.00 

0.6978 

0.04297 

0.6971 

3 

0.00 

0.6590 

0.02852 

0.6584 


extrapolation 0.00 0.7414 
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The values of the corresponding stress intensity factors are 
determined using the following relations: 

For Rice J integral: 


Kj 2 = 2 E ( Jj + yj Jj 2 - J 2 2 )/ {( 1+ v) (k+1 ) 

K n 2 = 4 E Jj / {(1 + V)(K+1) } -Kj 2 


3. 23 


For Bui J integral: 


K T 2 = 4 E Jj / {( 1 + v) ( k +1 )} 


K n 2 = 4 E J n / {(1 + v) (k+1 )} 


3. 24 


where V is the poisson's ratio and K is as defined earlier equals 
(3-4v) for plane strain condition. 

The corresponding values of the stress intensity factors 
from the J integral results are shown in Table 3.4, using 
relations 3.23 and 3.24, for Rice and Bui formulation, 
respectively. Also shown, the mode II stress intensity factor 
using the extrapolation technique. Both J integral formulations 
gave consistent results, which are on the average 9% below the 
extrapolation result. But looking at the results for different 
paths, the value of K changes by almost 5% between each path. 
The difference between the two J integral formulations is not yet 





obvious until the mi xed-mode loading condition is analyzed. 

For the mixed-mode loading condition, both cases analyzed 
above are combined together, see Fig. 3.14. The same dimensions 
and loading conditions are implimented in this example. Table 
3.5 summurizes the values of the J integrals from both 

formulations for this mixed-mode loading condition. Also shown 
are the corresponding stress intensity factors using Eq . 3.23 and 
Eq . 3.24 for Rice and Bui formulations, respectively, along the 

same three paths described above, see Fig. 3.11. While Bui’s 
formulation gives consistent results with the pure mode cases, 
Rice’s formulation gives completely different answers, especially 
for the Kjj stress intensity factor. A 33% difference in the 
value of the mode II stress intensity factor is observed, between 
the pure mode II and the mixed-mode case. This drop in the value 
of Kj is due to the incorrect integral formulation as 

discussed by Bui, Ref. [22]. Comparing the values of the mode II 
stress intensity factor between different paths, there is a 11% 
change in the J ^ between path 1 and path 3. The values of the 

stress intensity factors, and K^, from the extrapolation 

method are exactly equal to the pure mode results. 

Due to the above controversy, the extrapolation method is 
used in the application of the BIEM to mixed-mode failure of the 
inner raceway of the high speed bearings. 
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Table 3.5 Values of The J Integral Results with the 

Corresponding Stress Intensity Factors, for 
the Mixed-Mode Loading Condition, (?/t = 50). 


Rice Bui 


Path 

J i E 

j 2 e 

J 1 E 

J n E 


t 2 a 

2 

t a 

t 2 a 

+ *" a 


Eq . 

2.13 

Eq . 

2.16 

1 

4.1986 

1.8082 

3.7951 

0.6.35 

2 

4.2400 

1.7836 

3.7979 

0.6-22 

3 

4.1938 

1.6154 

3.7993 

0.3)o 


K I 

K II 

K I 

K. . 


t Ja 

t veT 

t/a 

tv a 


Eq. 3, 

.23 

Eq . 3, 

.24 

1 

2.0950 

0.4742 

2.0422 

0.6659 

2 

2. 1079 

0.4649 

2.0429 

0.6971 

3 

2.1049 

0.4217 

2.0433 

0.6584 

extrapolation 


2.0450 

0.7414 


CHAPTER FOUR 


APPLIED LOADINGS ON THE INNER RACEWAY OF BEARINGS 

As a first step towards attempting to determine the crack 
propagation rate and the time to failure of the high speed 
bearings, the loading acting on the inner raceway should be 
determined. The most important loadings are the tangential 
stresses due to the rotational speed and the press fit, and the 
next important loadings are the Hertzian stresses, which alter 
the tangential stresses with each passage of the heavily loaded 
roller. 

A.l Tangential Stresses 

The tangential stresses in the inner raceway are due to the 
rotational speed and the press fit of the inner raceway on to the 
shaft. Fig. A.l shows a simplified bearing configuration mounted 
on a hollow shaft. The solutions of the tangential stress in the 
inner raceway, modelled as a rotating disk is found to be, from 
Ref. [All: 

o ep = j(3+v)/8j P^ 2 (r^+r^+r^r^/r 2 - { ( l+3y)/(3+ • ) r 2 ) A.l 
where p is the density of the material, tu is the rotational speed 


87 


88 



Fig. 4.1 A Simplified Bearing Conf igurat icr. . 
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and r is the radial distance, and and are the inner and 
outer radii of the rotating disk. The most critical point for 
surface initiated crack propagation is at the outer radius of the 
raceway, r - r^. Substituting the value of for r in Eq . 4.1, 
the tangential stress at r = equals: 

° ee ■ |(3+v)/4 } pw 2 ^ 2 * + { (l-v)/(3+v)) r Q 2 * * * ] 4.2 

For typical aircraft engine bearings the ratio b-tween the outer 
radius to the inner radius of the inner raceway equals 1.13, Ref. 
[42]. Subtituting the value of r^ as a function of the inner 
radius of the raceway, r^ , and making use of the DN value, which 
is the product of the bearing bore in millimeter times the shaft 
rotating speed in RPM. The tangential stress at the outer radius 
will have the following relation assuming a poisson’s ratio of 
0.3: 


0 OQ = 3.32075 x 10 9 (DN) 2 (psi) 4.3 

00 

3 

where the density p Is assumed to equal 0.288 Ib.'in • A plot of 

the tangential stress versus the DN value Is shown In Fig. 4.2. 

Also shown Is an estimated critical crack length required for 

unstable crack growth. This estimate is derived for an edge 

crack in a semi-infinite plate with an applied normal stress 

equal to the tangential stress at the outer radius of the inner 


Ksi 
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DN, Million 


Fig. A. 2 Tangential Stress in the Inner Racewav t-.i Critical 
Crack Length versus DN Value in Millicr. 
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raceway. 

K t =1.12 o An / a tt A - A 
Ic tU! c 

where the value of K T Is the critical stress Intensity factor 
for unstable crack, growth. The value of used is equal to 
17.5 Ksi/Tn7 which is a typical critical stress intensity factor 
for M50 bearing steel. As seen from the equation, the critical 
crack length is inversly proportional to the fourth power of the 
DN value: 

a = 7.047 x 10 Zi4 (DN) _i+ (In) 4.5 
c 

The critical crack length required for unstable crack growth 
decreases rapidly for DN values greater than 1.5 trillion. For a 
DN value equals 3 million, the critical crack length is almost 
1/3 of the total raceway thickness. However, the effect on this 
critical crack length due to the press fit of the inner raceway 
onto the shaft ring should be estimated. 

The effect of the press fit is to increase the tangential 
stresses in the inner raceway causing the critical crack length 
to decrease. A typical interference fit of the inner raceway on 
the shaft is 0.00233 in (0.059 mm) at 0 RPM, Ref. (42], causing a 
uniform pressure p = 1,974 psi along the interference radius, for 
two cylinders having the following dimensions, r^ - -.6 in, r^ - 
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2.3 in and r' = 2.3023 in, r' = 2.0 in. The tangential stress 
o ’ i 

due to a uniform pressure distribution is given in Ref. [41]: 

°ee = p ( r i /(r o ' r i >) 1 1 + r o /r 1 4,6 

For r, = 2.3 in and r * 2.6 in the tangential stress is equal to 
i o 

16.178 ksi at r = r Q the outer radius of the inner raceway. But 
at increasing RPM the interference fit is reduced due to the 
radial displacement caused by the centripetal acceleration. 
Thus, a reduced pressure will develop between the inner raceway 
and the shaft which causes the interference fit to decrease with 
increasing speed: 

6 = 0.00233 - 1.0463 x 10 _16 (DN) 2 (in) 4.7 

new 

where $ is the new interference fit at a given DN value, 
new 

The pressure that develops from this new interference fit, will 
cause a tangential stress at the outer radius of the inner 
raceway which is given as: 

CJ 0e * 16, 178 - 7.265 x 10“ 10 (DN) 2 (psi) 4.8 

The net tangential stress at the outer radius of the inner 
raceway of the bearing, will be the sum of the two stresses due 
to the rotational speed and the press fit at a given DN value. 
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given as the sum of Eq . 4.3 and 4.8: 

°ee= 16, 178 + 2.59425 x 1 cf 9 (DN) 2 (psi ) 4.9 

and the critical crack length for the combined tangential 
stresses will be: 

a = 7. 77x 10 7 / (16, 178+2. 59425xl0" 9 DN 2 ) 2 4.10 

c 

A plot of the net tangential stress and the corresponding 
critical crack length acting at the outer radius of the inner 
raceway is shown in Fig. 4.3, The critical crack length 
calculated for a DN value of 3x10^ is much than the defect depths 
which are of the order of 0.004 in. Thus, the only cause for the 
crack to grow from a small surface furrow or a debris dent is 
through a fatigue process developed by the alternating Hertzian 
stress . 

4.2 Hertzian Pressure 

When two bodies are in contact, special attention must be 
paid to the deformation in the contact region. The original 
analysis of elastic contact stresses under static loadings was 
published by Hertz in 1882 [43]. In his honor, the stresses at 

mating surfaces of curved bodies are commonly called Hertzian 


C-^~" 



stresses 


• The actual load-carrying capacity of bodies in contact 
depends upon more than the elastic stress computed from Hertz 
equations* Surface sliding, lubricating fluid films and surface 
roughness should also be considered. In this section only the 
elastic solution under ideal surface conditions is assumed. 

Hertz's analysis showed that for two spherical bodies in 
contact the pressure distribution, which causes a uniform 
displacement field in the contact region, is represented by an 
elliptic pressure distribution, see Fig. 4.4, having the 
following equation: 



where p^ is the maximum pressure at the center of the contact 
region, and a and b denote the major and minor semiaxes of the 
contact surface for two spherical bodies in contact. The total 
load applied will be equal to: 


P = (2/3) Habp 0 


4. 12 


For two cylindrical bodies in contact, (roller bearings), the 
major semi-axis, a, goes to infinity and the pressure 
distribution will have the following form, see Fig. 4.3: 



4. 13 
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resulting in an applied total load which equals: 


P - (1/2) 77 b L Pq 


A. 14 


where b is the minor semi-axis and L is the length of the roller 
beari ng. 

To determine the dimension of the minor semi-axis, b, of the 
contact region, the curvature and the material properties of the 
two cylinders in contact should be taken into consideration. 
Following Ref. [44], b equals: 


b = 2 


p ( (l-v^) / Ej + d-v 2 ") / e 2 ) 


4.15 


TT L ( 1/rj + l/r 2 ) 


where r ^ , and are the radius of curvature, the poisson's 
ratio and the modulus of elasticity, respectively, for the two 
cylindrical bodies in contact. 

When the applied load for a given bearing configuration is 

known, the Hertzian stresses in the bodies can be easily 

determined analytically or numerically. Smith and Liu, Ref. [44] 

solved for the stresses for a semi -i nf i ni t e plate under a 

Hertzian pressure distribution using the integral equation 

technique. The stress distributions o , o , and c are given 

x y xy 

In Fig. 4.6, 4.7 and 4.8, respectively. Fig. *.7 shows the 
variation of the tangential stresses with roller position. As 
the roller gets closer to the line under considerations, the 
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stresses are highly compressive, and as the roller gets away the 
compressive stresses die out quickly. For a surface defect, this 
fast variations of the stress state superimposed on the 
tangential stresses due to the rotational speed and the press fit 
causes an alternating mode I stress intensity factor that can be 
a driving force for the fatigue crack growth. 

The second crack propagation mode is the shear mode, as seen 
in Fig. 4.8. The shear stress increases as the roller gets 
closer to the line under considerations, reaches its maximum at 
the edge of the loading and then decreases quickly. While the 
normal stresses decay relatively fast for distances of x greater 
than 5b, the shear stresses still have some effect at greater 
depth. 

Friction effect is sometimes superimposed on the normal 
Hertzian loadings for sliding contacts. The friction is usually 
represented as a fraction of the normal pressure acting in a 
direction opposite to the relative motion, given by Ref. [44]: 

q = f p 4.16 

where f is the coefficient of friction, and p in the normal 
Hertzian load. In general, in the absence of viscosity effects 
of the lubricant, friction changes little with speed. When 
viscosity effects comes into consideration, two types of behavior 
are observed; for mineral oils friction decrease with increasing 
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Depths x/b, when Friction is Present, (f=0.3) 
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speed, while for fatty acids the opposite is observed. At high 
speed, increase in friction is normally observed, Ref. (45]. The 
variations of the stress distributions under the Hertzian 
pressure with friction are shown in Fig. 4.9, 4.10, 4.11, for o , 
a and o respectively. The friction coefficient, f, used in 

y x y 

the calculation, equals 0.3. The net effect of the friction is 
to alter the symmetry of the normal Hertzian stress components; 
increasing the stress when the friction acts in the direction of 
the loadings and decreasing it for the other case. The most 
important change in the stress distribution is for o , which 
becomes positive when the loading is on the left of the line 
under consideration. For the shear stress, the negative 
component of the stress distribution increases, while the 
positive component decreases, Fig. 4.11. 

The next step in the analysis is to include the effect of 
the lubricating film since all high speed bearings are well 
lubricated. 

4. 3 Elastohydrodynami cs Lubrication 

For the completeness of the analysis the fluid film should 
be considered, in the presence of a surface defect or a surface 
furrow. 

The governing equation for the el astohyd rodynami c lubri- 
cation is based on Reynolds equation, which is a special case of 
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the Navler equation for fluid bodies. Assuming no side leakage, 
the two dimensional Reynolds equation is given by; 


d 

dy 


fph' 

tT 



dph 

12 v 

dy 


A. 17 


where p is the lubricant density 

n is the lubricant viscosity 
h is the film thickness 
p is the pressure 

v is the mean velocity of the lubricant (vj+V 2)/2 


The variation of the viscosity with pressure for liquids is: 


n = Hq * 




4.18 


where n .q Is the coefficient of the absolute viscosity at 
atmospheric pressure, and £ is the pressure-viscosity coefficient 
of the fluid. 

The density variation with pressure is: 


P = Pq [ 1 + ap/(l+bp) ] a. 19 

where Pq is the density at atmospheric pressure, and where a and 
b are pressure-density coefficients. 

The film thickness is given by: 
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h = h n + y^/2 (l/r.+l/r.) + u + d 4.20 

u 1 2 x 

where is the rainiinum film thickness, are the radii of 

curvatures of the two bodies in contact, u is the total elastic 

x 

deformation and d is the height of the surface defects in the 
contact region due to debris dents or surface furrows. 

The solution to Reynolds equation, Eq . 4.17, with the given 
viscosity and density as functions of the pressure can only be 
obtained numerically. The numerical technique used here is the 
combination of the finite difference method for the solution of 
Reynolds equation and the boundary integral method for the 
displacement calculations. 

Fig 4.12 shows the flow diagram for the numerical technique 

impl imented# An initial pressure profile is assumed. Then, the 

displacements for the given pressure profile are solved for, 

using BIEM program. The film thickness h, from Eq . 4.20, is 

calculated next. The new pressure profile is now determined 

using the finite difference method. Due to the high gradients of 
2 2 

dp/dy and d p/dy , a new smoother paramater, <p , is introduced, 
defined as: 

3 / 2 

4> = p h 4.21 

Note that p is small at large values of h and vice versa. This 
subsitution also has the advantage of eliminating all terms 
containing derivatives of h and p or h and <t> . Therefore, the 
Reynolds equation in terms of the new paramter c will be: 
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Fig. 4.12 Flow Diagram of the Elastohydrodvra-i; Lubrication 
Problem 
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Using standard finite difference representation, Eq. 4.22 can be 
rewritten as: 


a, <t> . + a- . 4> + a., . = b, 

1,1 i-1 2,i i 3,i l+l i 


4.23 


where 


1 » i 


3/ 

h 

i P 

a2 n 


i-h 


3, i 


h h 

r ' 


i+h 


2,i 


h 2 
= - -k < 

' P 

+ £ 

1 

y 

3 

— 4 

1 dh 

1 

1 §= 

► 

d 2 

l” 


i-h J 


2d 

[ n dy 

li + 1 n dy J 

i-l. 





< ph 

ph 



i+1 

i-1 

> 


d * nodal distance 


The system of equation of Eq . A. 23 is solved using the Gauss 
Siedel iteration method until $ does not vary. The relaxed Gauss 
Siedel method of iteration is given by: 


(I + w 1 aj) $ n+1 = |(l-w)I-wa2 1 a 3l^ >n + w a 2 ^ ^ 


where w is the relaxation coefficient. The Guass iteration will 
continue until <f)does not vary. 

The new pressure is calculated and compared with the 


no 


original pressure. If the difference is much bigger than a 
certain tolerance, a new displacement field will be calculated 
using the BIEM program. If not, the load generated from the 
pressure is measured and compared with the applied load required. 
If the two loads do not compare, a new film thickness h^ is 
estimated and the calculation is repeated until convergence. 

The boundary conditions for the Reynolds equation are: 

1) p = 0, for regions far from the contact area 

2) p is reset to zero at any position when the pressure is 

negative 

The second condition is commonly known as the Reynolds condition, 
which is based on obervations that a negative pressure cannot 
exist in the contact region. 

An initial run was performed using constant viscosity and 
density, with and without a surface furrow. The pressure 
profiles are shown in Fig. 4.13, for a defect depth equals 0.0030 
in. and a Hertzian load of 1500lbs. The pressure distribution 
drecreases sharply at the tip of the defect and increases to the 
original pressure distribution after the defect. As the depth of 
the defect increases the pressure in the pocket decreases to 
zero. Therefore, one can assume that with the presence of the 
surface crack the pressure of the lubricant in the crack can be 


taken as zero. 



Effect ot a bena 
Pressure Profile 
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4. 4 Initial Estimate of The Stress Intensity Factors 


As a first estimate, the stress intensity factors, and 
Kji, were evaluated using the superposi t i on method. This method 
consists of evaluating the stresses in an uncracked body and then 
applying the negative of the stresses for the uncracked body 
along the crack faces of the cracked body to determine the stress 
intensity factors. Fig. 4.14 represents a schematic diagram of 
the superposition method used. Fig. 4. 14. a shows a cracked plate 
with the tangential stresses and the Hertzian loading. This is 
equated to the sum of a plate with the same loading conditions 
but without a crack. Fig. 4.14.b. and a cracked plate loaded only 
along the crack faces with the negative stresses of the uncracked 
case. The stress intensity factors are then determined by 
summing the effect of every point load using available solutions. 
The stress intensity factors are given by: 


Kj. ' 2/ 

K n- 2/ 



(1 - (x/a) 


(1 - (x/a) 


2 


2 


) 


) 


F(x/a) dx 


F(x/a) dx 


4.25 


where F(x/a) is given in Fig. 4.15 from Ref. (-0], for a point 
load along the crack faces at distance x from the crack tip. The 
stresses used in this superposition method were the tangential 
stresses due to the press fit and the rotational speed combined 
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Fig. 4. 15 Solution for the st: 
a point load on the 
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with the Hertzian stresses for a given roller position and crack 
length. 

Fig. 4.16 represents the variation of the mode I stress 
intensity factor as a function of the roller position for 
different crack lengths, a/b. The value of the maximum Hertzian 
pressure, p^, used in this example, equals 180ksi, with a semi- 
major axis, b, equals 0.005 in. The change in the stress 
intensity factor decreases with increasing crack depth a. For 
small crack depths is negative, indicating crack closure, due 
to the high compressive Hertzian stresses. The variation of the 
mode II stress intensity factor with roller position is shown in 
Fig. 4.17 for different maximum Hertzian pressures Pq. The 
variation of the maximum stress intensity factor with crack 

depth for Pq equals 180ksi is shown in Fig. 4.18. The maximum 
mode II stress intensity factor increases sharply for small crack 
depths, reaches a maximum value, and then decreases slowly as the 
crack depth increases. These results show the variation of the 
stress intensity factors with each passage of a roller. As seen 
from the figure, the stress intensity factors are well below the 
critical stress intensity factor, K^, of M50 bearing steel, for 
a crack depth less than 0.06 in. Therefore, the only criteria 
for the failure of the inner ring from a very small surface 
defect is through a fatigue process due to the combined 
alternating stress intensity factors, A and as s h own 


in the next chapter 


116 


Although this method gives accurate solutions for the stress 
Intensity factors for straight cracks, a more generalized method 
should be used to determine the complete stress and displacement 
distributions around any arbitrary branched crack. The 
generalized method that will be used to determine the stress and 
displacement fields as well as the stress intensity factors, for 
any arbitrary oriented crack, is the BIEM as will be seen in the 
next chapter. 
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Fig. 4.18 Maximum YL Stress Intensity Factor as a Function of 
Crack Length 



CHAPTER FIVE 


MIXED-MODE FAILURE ANALYSIS OF THE INNER RACEWAY 
OF HIGH SPEED BEARINGS 

Chapter One described the different failures observed in the 
inner raceway of bearings, and especially the new catastrophic 
failure recently detected in test rigs for high speed aircraft 
engine bearings. This chapter utilizes the mixed-mode failure 
criteria of chapter Two to analyze the causes that lead to the 
brittle failure of high speed bearings, using the boundary 
integral method described in chapter Three and the major loading 
conditions summarized in chapter Four. The interaction of the 
Hertzian stresses and the tangential stresses due to the 
rotational speed and press fit, in the presence of a surface 
defect are first described in terms of the stress intensity 
factors, Kj. and . Then, the crack propagation direction is 

determined using the different crack growth criteria described in 
section 2.4. Finally, the time to failure is estimated using 
fatigue failure data of M50 bearing steel. 

5.1 Stress Intensity Factors 

The stress intensity factors are determined for a typical 
roller bearing used as a support for the main shaft of aircraft 
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Table 5*1 Typical Roller Bearings for Turbine 
Engine Mai n Shaft . 

-Roller Diameter 0.5 in. 

-Roller Length 0.57322 in. 

-No. of Rollers 28 

-Shaft Inner Radius 2.0 in. 

Outer Radius 2.30233 in. 

-Inner Raceway Inner Radius 

Outer Radius 

^ T> T ^ 

Outer Radius 
-Shaft Speed 25,500 RPM 
-Load on most heavily loaded roller 500 lbs., 1,000 lbs., 

1,500 lbs., 2,000 lbs. 

-Interference fit of inner raceway 0.00233 in. 


2.300 in. 
2.600 in. 

3.35 in. 


onto the shaft at 0 RPM 
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engines. The dimensions of the bearing and the shaft are shown 
in Table 5.1. A two dimensional analysis is performed assuming a 
plane strain condition, representing the mid-section of the inner 
raceway. This approximation is assumed sufficiently good to 
identify some mechanisms that lead to the catastrophic failure of 
the inner raceway. The inner raceway was modelled as a plate 
with width equal to the inner raceway thickness and length twice 
the width, corresponding to the distance between two rollers. 
The curvature of the raceway was ignored since the ratio of the 
raceway thickness to the outer radius of the inner raceway is 
smal 1 . 

The tangential stresses applied at the edge of the plate are 
calculated from the solution of a rotating disk at 25,500 RPM 
superimposed onto the tangential stresses of the press fit at the 
same rotational speed. The calculated tangential stress at the 
middle of the inner raceway was 42,300PSI. The variation of the 
tangential stress along the radius was approximated by a linear 
curve, decreasing to about 7% of the mid value. 

A reaction pressure load was applied at the bottom of the 
plate to simulate the contact stresses transfered from the inner 
raceway to the shaft. The corresponding load of this pressure 
will be equated to the applied Hertzian load. Different Hertzian 
loading condition are assumed giving rise to different maximum 
pressures p^ and contact lengths b. 

Table 5.2 summarizes the four Hertzian loading conditions 
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Table 5.2 The Geometry of the Hertzian Pressure 
distribution for different loads. 


P lbs. 

b in* 

Pq ksi 

500 

0.00392 

142. 

1,000 

0.00554 

200. 

1,500 

0.00679 

245. 

2,000 

0.00780 

283. 
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used in the analysis- Fig. 5-1 shows the geometric modelling of 
the inner raceway and the number of nodes used in each 
subdivision for BIEM analysis. 

The stress intensity factors are determined using the 

extrapolation method described in section 3.5. The shear 

stresses o and the normal stresses 0 are multiplied by the 
xy y 

square root of the distance to the crack tip and plotted versus 
that distance- Their intercepts with the axis at r equals zero, 
are proportional to the stress intensity factors and 

respectively. The stresses used in the extrapolation method were 
taken between 0.2a and 0.5a ahead of the crack tip. 

To simulate the passage of each roller, the Hertzian loading 
is moved incrementally along the top face of the plate. The 
dynamic effect was ignored assuming steady state solution. 

A typical variation of the mode I stress intensity factor 
with roller position and crack depth is shown in Fig. 5.2, for a 
Hertzian load equals 1,500 lbs. As seen from the figure, the 
compressive stresses of the Hertzian loading cause the decrease 
of the constant value of the mode I stress intensity factor due 
to the tangential stresses of the rotational speed and the press 
fit. For very small cracks, is constant and then decreases 
rapidly as the roller approaches the crack and increases rapidly 
again after the roller moved away from the crack. The negative 
Kj values observed indicates crack closure. For large cracks the 
rapid decrease in diminishes and practically disappears since 


(in. 
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Fig. 5.2 Variation of K Stress Intensity Factor with 

Roller Position as a Function of Crack Length, 
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the crack tip at these lengths has passed outside the highly 
compressive Hertzian stress field* 

These variations in can result in crack growth for small 
cracks. However, there is another stress intensity factor namely 
the shear mode or mode II, which may play a more important role 
in crack propagation. The variation of the mode II stress 
intensity factor Is shown in Fig. 5.3 for the same Hertzian 
loading condition as above. The value of starts at zero when 
the roller Is far from the crack tip. As the roller approaches 
the crack, K starts to decrease until it reaches a minimum, as 
the Hertzian nressure Hi strihnt 1 on rearhea the rrark. then 
increases rapidly back to zero. When the roller starts to move 
to the other side of the crack it changes sign and starts to 
increase to the maximum value and then decreases back to zero. 
For small cracks the change in is very abrupt, but for deeper 
cracks, the change in Is more gradual. The change in 

persists at greater depths than K^. 

The variations of the stress Intensity factors with roller 
position for different Hertzian loading are shown in Fig. 5.4 and 
5.5 for Kj and K^, respectively, for a crack length equals 
0.008in. As seen in the plots, the greater the Hertzian load the 
greater the change in the stress intensity factors. 

Fig. 5.6 represents the maximum absolute value of for 
different applied Hertzian loads. The values of increase 
until reaching a maximum at depths which correspond to the 
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Hertzian Load 
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Fig. 5.4 Variation of Kj Stress Intensity Factor with 
Roller Position as a Function ol the Hertzian 
loads, (a=0 . 008 in . ) 


Hertzian Load 


130 



h C 


AJ 

c 

•M 

<0 

5 

•H 


N 

Ui 

AJ 

O 

1* 

AJ 

OJ 

u 

X 

to 


u. 

<u 


X 


AJ 

AJ 



U-l 

03 

o 

c 



c 

i_J 

o 

c 


1— 1 

AJ 


u 

CD 

c 

03 

3 

<u 

&u 

U 


AJ 

co 

C/3 

03 


*— i to ^ 

w c 
c 

o oo 
u-i ' r " t O 
0^0 


M II 
O to 
a* w 

v- 
<u 


V- 

to 

> 


m 

in 

U* 


( -UT/ -TS^) 


Us 


Loads 


Hertzian Load 


131 



( *_uy ‘ ts^) 


Maximum Absolute Value of K S:ress Intensi 
Factor for Different Hertzian L^ads as a 
Function of Crack Length 
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maximum shear stress and then decrease slowly. This is due to 
the fact that although the shear stress decreases as we move down 
from the surface causing to decrease, the increase in crack 

length causes to increase. 

When friction is added the stress intensity factors are 
modified as shown by the stress distributions given in section 
4.4. The variation of the stress intensity factors and 

for different friction factors, f , ( f = 0.0, 0.05, 0.1, 0.2, 

0.3), are plotted in Fig. 5.7 and 5.8, for a crack length equals 
0.002in. As the roller approaches the crack the friction forces 
act in a direction to close the crack, causing to decrease 
even more, but as soon as the roller passes the crack, the 
friction forces act now to open the crack, causing an increase in 
Kj. As the roller moves further out, the effect of friction 
disappears. As for mode II, as shown in Fig. 5.8, increases 

in absolute value as the roller approaches the crack, and then 
decreases as the roller moves to the other side of the crack 
since at that point the shear stresses are positive while the 
friction forces are still causing a negative shear. 

Note that the stress intensity factors for the crack length 
considered are below the critical value, K^, anc * w ^ en * s 

maximum is zero and when ^ is maximum, is almost 

minimum. This lower than critical stress intensity factor and 
phase shift causes a complex fatigue mixed-mode loading condition 
that can only be analysed using the crack propagation criteria 
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of Roller Position for Different Sliding Friction 
Factors, f, (a-0.002in.) 
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discussed in chapter TWo. 

5*2 Crack Propagation Directions 
a) Static Crack Growth 

To determine the crack propagation direction, an effective 
stress intensity factor should be considered to combine the 
mixed-mode loading effect observed with each passage of the 
roller. The criterion, that will be mainly emphasized, is the 

1 * * * • i ^ _ _ 1 . . _ . _ »» 

O fV CALCIlOXUll 1 ^ L L C L L i L C L X W II , Q 1 O \J l\ uu w u Cl J l • i u u i g j 

Release Rate", as discussed in section 2.2. Therefore, the 
stress intensity factors for every roller position, will be used 
to determine the crack extension forces due to the tangential 

stresses, G , and due to shear, G , for every angle 0 around 

ou r 0 

the crack tip, using Eq • 2,10. 

A typical variation of the maximum crack extension forces, 

G ^ and G^, with roller position, are shown in Fig, 5.9 and Fig. 
5.10, respectively, for a crack length equals O.OQ^in. The crack 
extension forces due to the tangential stresses have a constant 
value when the roller is far from the crack, (y^/b is large). 

G^ starts to decrease with the approach of the highly 
compressive Hertzian stresses, then starts to increase to a new 
maximum due to the sharp increase in the mode II stress intensity 
factors present when the Hertzian pressure is at the edge of the 
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0 .004 in. 
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crack, (y 0 /b=1.0). When the roller is positioned on top of the 
crack tip, (y 0 /b=0.0), the shear effect goes to zero, and the 
crack is completely closed, causing the crack extension forces to 
be zero, assuming no lubricant had leacked in the crack. The 
same behavior is observed for the crack extension forces due to 
shear, as in Fig. 5.10. 

With the presence of the Hertzian loading, two competing 
mechanism exist; 1) due to the compressive Hertzian stresses, the 
crack extension forces are decreased, and 2) due to the shear 
stresses, the crack extension forces are increased. 

The angle along which the maximum tangential crack extension 
force acts, as a function of the roller position is shown in Fig. 
5.11 and 5.12, for the maximum direction of G and G , 

ee re 

respectively. The angle starts at zero for the direction along 
which the maximum tangential crack extension forces acts, 
implying is dominant, then starts to increase as the roller 
approaches, until reaching a maximum value of 70.53 degrees, 
indicating a pure shear mode, at y^/b^l.C. The maximum shear 
crack extension force acts in the opposite direction; acts along 
70 degrees for the pure mode I case, for large >v'b, and at zero 
degrees for the pure mode II case, at yp/b=0.0. 

Fig. 5.13 shows a comparison of the values of the different 
mixed-mode crack propagation criteria; the maximum tangential 
stress, the minimum strain energy density factor, and the J 
integral approach as defined in section 2. A. The tangential 
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Fig. 5.11 Variation of Direction, 9, of the Maximum Tangential 
Crack Extension Forces as a Function of Roller 
Position, ( a=0 . 004in . , p=1500.1ts.) 
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FiR. 5.12 Variation of the Direction, 0, ofthe Maximum Shear 
Crack Extension Forces as a Function of Roller 
Position, (a=0 .004in. , P=1500. lbs.) 
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Fig. 5.13 Comparison between Crack Growth Criteria as a Function 
of Roller Position, (a=0.004in. : P=1500. lbs.) 
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Fig. 5.14 Comparison between the Orientation of Crack Growth for 

Different Crack Extension Criteria, (a=0.004in., P-1500, lbs. 
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stress criterion corresponds exactly to the tangential crack 

extension forces, G . All the criteria show the same trend of 

0 6 

decreasing in value as the Hertzian stress approaches and 
increase sharpely due to the shear effect. The main difference 
lies in the angle along which the maximum values act, as seen in 
Fig. 5.14. The minimum strain energy criterion shows a maximum 
at 83.61 degrees while the J integral at 0 degree, and the 
tangential stress at an angle of 70.53 degrees, at y Q /b=1.0. 

Note that the maximum mode II stress intensity factor occurs 
with negative K values for small cracks, indicating crack 
closure. It had been observed that when mode II is associated 
with compressive normal stresses, friction occurs along the 
contact crack faces, thus reducing the effective mode II stress 
Intensity factors, Ref. [47]. Swedlow, Ref. [48], proposed a 
simple coulomb friction along the contacting faces of the crack, 
acting in the opposite direction of the shear and proportional to 
the normal compressive stresses. Since the crack modelled In the 
BIEM mesh has a finite width, the compressive normal stresses on 
the crack faces are absent. As a result, the effective mode II 
stress intensity factor, was approximated by: 


K tt „ = K TT + f K t if L < 0 and K TT i 0 
lief f III I II 


5. 1 


K IIeff ~ K II 


If kj > 0 
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Fie. 5.15 Variation of the Maximum K n for Different Friction values 
along the Crack Faces as a Function of Crack Length. 
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where f takes the opposite sign of K^, since the friction acts 
to the decrease the shear. Although this model is crude, it is a 
good first approximation of the friction effect. The effect of 
the friction model used on the maximum mode II stress intensity 
factor is shown in Fig. 5.15 as a function of the crack length 
for different friction coefficients. ^ decreases with 

increasing friction factor until the crack length reaches 0.012in 
where K^. becomes positive, and no variation in is observed. 

The effect of the sliding friction, (friction between the 
roller and the inner raceway), on the crack extension forces is 
shown in Fig. 5.16 and 5.17. Fig. 5.16 represents the variation 
of the tangential crack extension forces as function of the 
roller position, for a crack length equals 0.03-in. For y^/b 
less than zero the crack extension forces increase with 
increasing frictions, while decreasing with increasing friction 
for yg/b greater than zero. This behavior is due to the friction 
forces which act to open the crack on one side of the crack and 
to close it when it Is on the other side. The crack extension 
forces due to shear behave in almost the same manner, as seen in 
Fig. 5.18. 

The third mechanism that also alter the variation of the 
crack extension forces, is the lubrication effect when the roller 
is on top of the crack, causing a normal stress distribution 
along the crack faces. This mechanism can modify the crack 
extension forces if the lubricant can penetrate into the crack 
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for Different Sliding, Friction Factors, as a 
Function of Roller Position, (a=0.004in., P=1500. lbs. 
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when the roller is on top of the crack tip. This mechanism 
involves a complex fluid-solid interaction problem, which has to 
include the viscous-drag and the surface tension forces to 
determine how far the lubricant can penetrate into the crack. 

But as the roller approaches the crack, the normal 
compressive stresses, o^, try to close the crack and push the 
lubricant out, before the roller gets on top of the crack. For 
longer cracks, the commpressive stresses are limited to a small 
distance close to the surface, and at greater depths the 
stresses, o , are in tension, due to the tensile stresses of the 
rotational speed. This effect causes the lubricant to be caught 
inside the crack and to pressurize it. 

A uniform pressure was applied along the crack faces, when 

the roller is on top of the crack line, (yQ/b*0.0). The pressure 

assumed equals the maximum Hertzian pressure, p^. The variation 

of the mode I stress intensity factor as function of crack depth 

is shown in Fig. 5.18, for y Q /b=0.0. The values of are now 

positive, due to the high pressure applied along the crack faces, 

as compared to the negative values when no lubricant is assumed 

leaked inside the crack. The values of increase rapidly with 

increasing crack length, and reache the critical stress intensity 

factor K_ at a crack depth equals 0.007in. This model is not 
Ic 

realistic since the crack is closed as the roller approaches the 
crack line. If the crack is filled with lubricant, the high oy 
stresses will push the fluid outside. So assuming that there is 
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Fig. 5.18 Variation of Kj Stress Intensity Factor as a 

Function of Crack Length when Pressure is Applied 
Along the Crack Faces, (ypVb^O.O) 
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no lubricant leakage, Is a good approximation for very thin 
cracks. 

Note that the results for the crack extension forces and the 
crack growth directions are for static loadings at one particular 
roller position. The crack extension forces given above are well 
below the critical value for unstable crack growth, for crack 
length below 0.05 in. assuming no-lubricant leakage. Thus, a 
fatigue process has to take place for a small crack to grow from 
the surface to the critical crack length. 


b) Fatigue Crack Growth 


The fatigue process is governed by the maximum changes in 
the crack extension forces, AG^, and the mean value G^, as seen 
in section 2.5. For a given crack length, a, the maximum and 
minimum values of the crack extension forces are monitored for 
every angle 0, around the crack tip, as the roller approaches the 
crack and moves away from it. The direction, 0 , which has 

the maximum change in the crack extension force, will be the 
direction of the crack growth. This model is based on the 
assumption that an element oriented in the direction of the 
maximum change in accumulates the largest amount of damage 

with each passage of the roller, making this element the first to 


open up and join the main crack. 

Fig. 5.19 and 5.20 show the values of the maximum changes in 
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the crack extension forces, AG and AG , as functions of crack 

re ee 

depth. The maximum change in the tangential crack extension 
forces keeps on increasing with increasing crack length, 
Fig. 5. 19, while the change in the shear crack extension force 
increases sharply until reaching a maximum value then decreases, 
gradually, Fig. 5.20. 

The actual directions of the crack growth could be plus or 
minus the presented values since we have complete reversal of the 
stress intensity factors when the roller cross from one side of 
the crack to the other. This symmetry will break down only when 
sliding friction is present. 

The angle along which the maximum AGqq acts, starts at a 
constant value of 70.53 degrees for very small cracks and high 
Hertzian loads, and then starts to decrease gradually with 
increasing crack length, Fig. 5.21. This constant value of the 
direction of the maximum AGqq is due to the pure mode II stress 
intensity factor associated with negative K^. For very small 
cracks and low Hertzian loadings, the change in is negligible 

compared to the change in Kj, due to the high compressive 
Hertzian stresses and low shear stresses, causing the maximum 
change to occur along zero degrees (as seen in Fig. 5.21 with a 
Hertzian load of 5001b. and a crack length less than 0.004 in.). 

The angle along which the maximum change in the shear crack 
extension forces acts, starts at zero and increases gradually 
with increasing crack lengths, Fig. 5.22. The crack length where 
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Fig. 5.19 Variation of the Maximum Change in the Tangential 

Extension Forces as a Function of Crack Length, for 
Different Hertzian Loads. 
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the angle starts to deviate from zero Increase with the increase 
In the Hertzian load. 

To determine which mode is the crack driving force, (shear 

growth or tensile growth), the values of and A G are 

Uo r u 

compared with the threshold values for the material under 

consideration. Fig, 5.23 gives the regions for different crack 

growth modes. For values of AGg 0 greater than 0.49 psi-in, Ref. 

[49], the tensile mode is dominant and will control the crack 

growth direction and rate. For value of AGggless than AG^^, 

two regions are considered depending on the value of * For 

r 0 

values of AG . less than AG there Is no crack growth, while 

r0 r0th 

for values above AG . , , the shear mode Is the crack driving 

r 0 1 n 

force. 

Taking into consideration the above crack growth regimes, it 
is obvious that for crack depths less than 0.002in, the crack 
growth is governed by shear, while for a crack depth equal and 
above 0.002in the tensile crack growth is dominant. 

When sliding friction is present, (friction between the 
roller and the inner raceway), the change in the crack extension 
forces increases with the increase of the friction factor. The 
main effect of friction is to make the crack extension forces 
non-symmetric, causing the crack to grow mainly in the direction 
opposite to the rolling direction, as observed experimentally, 
Ref. [50]. Table 5.3 presents the variation of the crack 
extension forces and the directions they act on and a Hertzian 








158 


Table 5.3 Variation of the Maximum Changes in Gqq , 
G and their corresponding angles, for Different 
Friction Coefficient for two crack Lengths, 

P = 1,500 lbs. 


a a 


f 


0.002 

in. 



0.004 

in. 


AG ee 

0 

AG r0 

6 


0 

hC r6 

0 

0.00 

0.81528 

71. 

0.61152 

0. 

2.2626 

71. 

1.69716 

0. 

0.05 

0.97454 

71. 

0.73098 

0. 

2.5972 

71. 

1.94813 

0. 

.1 

1.1424 

71. 

0.85689 

0. 

2.9604 

71. 

2.22055 

0. 

.2 

1.5231 

71. 

1.14247 

0. 

3.7565 

71. 

2.81766 

0. 

.3 

1.9567 

71. 

1.46769 

0. 

4.6473 

71. 

3.48578 

0. 
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Table 5.4 Variations of the Maximum Change in Gqq 

and G Assuming a Coulomb Friction Along the 

ft) 



Crack 

Faces, a * 

0.004 in. and 

P = 

f 

AG ee 

e 

AG 

re 

e 

0.000 

2.2626 

71. 

1.69716 

0 . 

0.005 

2.2484 

71. 

1.68649 

0 . 

0.010 

2.2342 

71. 

1.67585 

0 . 

0.030 

2.1780 

71. 

1.63365 

0 . 

0.050 

2.1224 

71. 

1.59198 

0 . 

0.100 

1.9867 

71. 

1.49016 

0 . 

n onn 

1 7?«6 

71 

1 7Q661 

n_ 

0.300 

1.4885 

71. 

1.11652 

0 . 

0.500 

1.0622 

71. 

0.79672 

0 . 

0.800 

0.6758 

0 . 

0.41796 

0 . 
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loading equal to 1,500 lbs., for different friction coefficients. 
The effect of friction Is to Increase the maximum change of the 
crack extension forces with the Increase In the friction value. 
The direction, where the maximum changes in the crack extension 
forces occurs, does not change with friction. 

Table 5.4 shows the variations of the maximum change, of the 
crack extension forces for a crack length of 0.004ln. when 
friction along the crack faces are considered, as modelled in Eq. 
5.1. The effect of friction is to reduce the mode II stress 
intensity factor when the crack is closed. The crack extension 
forces decrease with increasing friction, but the crack growth 
direction does not change, until large friction is assumed, 
causing the shear to stop. This vanishing effect of the shear 
mode is observed at a friction factor of 0.8 for a Hertzian load 
of l,500lbs, and a crack length equal to 0.004in. 

c) Crack Kinking 

As seen from the above calculations, the crack tries to 
change the direction of growth. To simulate this effect, an 
original crack, a^ of 0.004in., is extended in the direction of 
the maximum change of the crack extension forces for three 
different new crack ligament lengths, a 2 > as seen in Fig. 5.24. 
The new stress intensity variations are determined along the new 
crack ligament, for different roller positions. 
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Table 5.5 Variation of the Maximum Crack Extension 
Forces for Different Hertzian Load and 
Kinked Ligament Length. 


P lbs. 

0 2 

in 

a 2 in 

> 

O 

CD 

CD 

6 3 




^0.001 

0.61261 

20 

500 

69 

0.004 < 

0.002 

0.65259 

18 




. 0 . 003 

0.94323 

-3 




j 0.001 

0.57596 

21 

1500 

71 

0.004 <0.002 

0.62055 

19 



l 0.003 

1.20560 

0 
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and the Angle they act on as a Function of Roller Position, 
(a *0.004in. , a =0.002in., 6 =71°, P=1500 lbs.) 
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Table 5.5 summarizes the loading conditions and crack 

geometries for the kinked crack for different Hertzian loads. 
The new crack Increment was extended In the direction of the 
maximum change In the tangential crack extension forces, AG^q . 
Three different crack increments were used,( i.e. a 2 “ 0.001 in., 
0.002 in. and 0.003 in.), to determine if the crack increment 
length will affect the direction of crack propagation. Figs. 
5.25 and 5.26 show the variations of the crack extension forces 
as a function of roller position for a Hertzian load of 1,500 

lbs. and a kinked ligament of 0.002in., for G and G respec- 

00 r0 

t" T7 O 1 \7 . UKon f ho r n 1 1 or “to far f r Am h Ko H no t-Ko 

tangential crack extension force starts out around 0.6 psi-in., 
and decreases to zero when y^/ b reaches -2. Due to the high 

shear stresses, the value of Gqq increases sharply to 0.24 

psi-in. and then decreases to 0.11 psi-in and then increases back 
up to a value of 0.6 psi-in. 

The direction of the maximum Gqq starts at 19 degrees then 
decreases to 0 degrees. As soon as the roller is on the other 
side, yg/b is greater than 0., the maximum Gqq occurs almost in 
the original direction then starts to decrease back to 19 
degrees. 

The maximum changes in the crack extension forces and the 
angles they act on are shown in Table 5.5. All the changes of 
G 0 q are above ^gg t h of 0*49 psi-in. for M50 steel, indicating 
aitensile mode growth. As seen from the results the crack tries 
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Fig. 5.27 Experimental Crack Profile, Ref. [51] 
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Table 5.6 Variation of the Direction of Crack 
Growth with Every Crack Extension 
Ligament, P=l,5001bs. 


i 

V 

AG ee 

9 i+l 

1 

0.004 

2.2626 

71. 

2 

0.002 

0.62055 

19. 

3 

0.002 

9.1371 

65. 

4 

0.002 

2.2293 

14. 

5 

0.002 

8.8994 

45. 
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Fig. 5.28 Simulated Crack Growth Profile for a 
Roller Bearing with a Hertzian equal 
to 1500. lbs. 
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to realign itself in the original direction, going back halfway. 
If the crack is extend further (a 2 - 0.003in.), the angle of 
crack propagation return to zero. We can state that: as soon as 
the crack grows from the vertical direction along the maximum 

angle of AG , the loading cycle changes causing the crack to 

0 0 

turn back to its original direction. This may explain the step 
like crack growth observed in the experiment of Bamberger at DN 
value of three million, see Fig* 5.26 from Ref. [51]. 

To simulate the experimental crack profile, the crack 

geometry was extended along the new maximum change of G , for 

00 

a constant ligament length, a^ equal to 0.002 in. This constant 
ligament length is arbitrary. 

Table 5.6 shows the values of the maximum change in G and 

00 

the angles where they act on, measured from the vertical, 
assuming equal crack increments of 0.002 in. Fig. 5.27 is a plot 
of the crack profile from the boundary integral simulation. One 
can observe some similarity between the actual crack profile of 
Fig. 5.26 and the simulation in Fig. 5.27. Qualitatively, one 
can observe the step-like variation of crack orientation with 
crack depth. 

The next step is to determine the time to failure of a small 
induced crack on the surface of the inner raceway of a high speed 
bearing, as will be discussed in the next section. 
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5.3 Time to Failure 


The time to failure of the inner raceway is estimated using 
the experimental results on the crack propagation rate for M50 
steel. Fig. 5.29 shows the crack growth rate as a function of 
the stress intensity factor range. The slope equals approxiraatly 
3. The shape of the curve is similar to the last segment of the 
da/dN verus AK^ curve, or the last stage in the crack growth 
rate, (Fig. 1.4). 

The time is estimated using Eq . 2.27 for tensile mode 

growth. The different contants, n^ and c, for different values 
of n£ are shown in Table 5.7, since no exact value of n^ is 
available. The life is estimated at 25,500 RPM starting with an 
initial crack of 0.002 in. N is determined by integrating Eq . 
2.27 with respect to da: 


N = 



da 


f(AG 00 , R) 


5.2 


The function f(AG ^ R) is approximated using least square of the 
values of AG as function of crack length a from Fig. 5.19. 

OO 

The variation of the crack length as function of the number 
of cycles is shown in Fig. 5.30. The main characteristic of the 
curve is the slow propagation phase for small crack length 
followed by rapid growth to failure. 
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Fig. 5.29 Crack Propagation Rate as a Function of 
the Stress Intensity Factor Range for M50 
Steel, Ref. [49] 
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Table 5.7 Constant used In Eq . 2.27 for M50 Steel 


case 

n l 

n 2 

C 

1 

1.433 

0.0 

1.5x10 

-7 

2 

0.82 

0.5 

8.4x10 
„ , -6 

3 

0.76 

1.0 

2.3x10 


173 


The variations of the number of cycles necessary to 
propagate a crack to failure as a function of the applied 
Hertzian load are shdwn in Table 5.8 with and without friction 
along the crack faces. The effect of friction along the crack 
faces is to increase the life to failure due to the decrease of 
the change in the crack extension forces with an increase in the 
friction factor. Fig. 5.31 shows the Hertzian load versus number 
of cycles to failure for an .original crack length of 0.002 in. 
The behavior of the curve, which shows an increase in life for 
decrease the Hertzian load on bearings, is similar to 
standardized a-w curves. 

A direct comparison with the experimental fatigue life of 
bearings at DN value of three million was not possible since the 
crack growth rate data available are for M50 steel tempered once 
while the bearings actualy tested were tempered five times. It 
is well known that the crack growth rate changes with the number 
of temperings. The second variable that affected the comparison 
is the lack of crack growth rate data of the shear crack growth 
mode. For very small cracks, the crack driving force, as seen 
earlier, is the the shear crack extension force, A q 

r e 

Nevertheless, one can observe the very fast time to failure of 
the inner raceway which is in the order of 10 minutes for a 
original crack length of 0.002 in. with a Hertzian load of 1500 
lbs. and a rotational speed of 25,500 RPM. This estimate is 
found assuming that the maximum loaded roller passes over the 


Crack Length a (in. 
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Table 5.8 Number of Cycle to Failure for Different 
Hertzian loads and friction values, f, 
along the crack faces. 


P lbs. 


500 1000 1500 2000 


case 

1 

f=0.0 2 

3 


0.946x10'* 

0.796xl0 5 

0.991xl0 5 


0.500x10"* 

0.475xl0 5 

0.435xl0 5 


0.296xl0 5 

0.287x105 

0.215xl0 5 


0.199xl0 5 

0.176xl0 5 

0.098xl0 5 


1 

f=0 .05 2 
3 


0.509xl0 5 

0.482xl0 5 

0.443xl0 5 


0.306x10"* 

0.295xl0 5 

0.226xl0 5 


0.208xl0 5 

0.177xl0 5 

0.107xl0 5 


1 

f=0. 1 2 

3 


0.518x10"* 

0.490xl0 5 

0.451xl0 5 


0.318xl0 5 

0.304xl0 5 

0.235xl0 3 


0.220xl0 5 

0.187xl0 5 

0.1l7xl0 5 



Hertzian Load P, (lbs. 
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crack only once every cycle. This results can explain the fast 
crack propagation of an induced surface crack on the inner 
raceway of bearing with a DN value of three million. 

5.4 Conclusion 

The main conclusion to be drawn from the preceeding analysis 
is that the alternating mixed-mode loading is the driving force 
behind the fast growing crack in high speed bearings. This 
alternating mixed-mode is due to the superposition of the 

an c t" t* o c c f -i o 1 rl rtnfn tho t* a ncron t* { a 1 cfrpCQoc rluo fn fhp 

rotation and press fit. The crack growth direction is shown to 
propagate in a step like fashion alternating between the vertical 
and an angle approaching 80 degrees as observed experimentally in 
Ref. [51]. 

The maximum change of the crack extension force components, 

AGqq and AG^are good measures of the crack propagation driving 

forces and direction of crack growth. 

The time to failure of a surface defect was determined to be 

very short, in the order of 10 minutes, which did not compare 

well with the experimental values of 2 hours for ball bearings. 

Crack growth rate in the shear mode growth region, ( region 

where AG is dominant) should be available to better predict the 
r“ 

fatigue life of the inner raceway and any other component under a 
general mixed-mode loading condition. 


178 

It should be emphasized that many simplifications were made 
in the present analysis. In particular the three-dimensional 
effects, thermal effects, residual stresses, dynamics, damage 
accumulation ahead of the crack tip, and finally the randomness 
inherent in the system, were not included. Although the results 
may not be the final answer to the fast failure of the inner 
raceway of high speed bearings at DN values above three million, 
some more understanding on the mechanism that leads to this new 
catastrophic failure was obtained. 


CHAPTER SIX 


CONCLUSIONS 


6. 1 Summary 


The goal of this study was to analyze the causes that lead 
to the brittle failure of the Inner raceway of high speed 
bearings. The analysis has shown the following conclusions: 

- The critical crack length required for unstable crack 
growth is inversely proportional to the fourth power of DN. 

- A mixed-mode fatigue crack growth is the cause that leads 
a small surface defect to propagate inward until reaching the 
critical crack length. 

- The mixed-mode loading is associated with the interaction 
of the compressive normal stresses and the shear stresses with 
the passage of the Hertzian pressure, and is superimposed on the 
tensile stresses of the rotational speed and the press fit of the 
inner raceway onto the shaft. 

- The maximum changes in the crack extension forces (AG , 

0 0 

AG ) are good measures of the mixed-mode fatigue driving forces. 

^ 0 

- Improvement in the fatigue properties of the bearing 
materials are needed to increase the threshold values of the 
crack extension forces in shear as well as in tension, (increase 


the values of 


AG , , and 
r 0th 


AG eeth ) * 
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- For very small cracks the shear crack growth Is the 
dominant crack extension driving force ( AG^). 

- The crack was shown to grow In a step like fashion as 
observed experimentally. 

- The effect of sliding friction is to make the crack grow 
in the direction opposite to the rolling direction. 

- The effect of friction along the crack faces is to 
decrease the shear mode stress intensity factor in the presence 
of high compressive normal stresses. 

- The time to failure was estimated to be very short in the 
order of 10 minutes. 

6.2 Further work 

- A three dimensional analysis for roller bearing as well as 
ball bearings should be considered. 

- The effect of the lubricant should be considered in 
greater details with viscosity pressure and temperature 
variations. 

- Residual stresses and material nonlinearity can be the 
cause of crack arrest and slower crack propagation rate as 
observed experimentally rather than the very fast growth measured 
from the elastic analysis. 

- More understanding of the crack growth behavior in shear 

growth mode region, where AG is dominant, is required. 

r 0 
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- Dynamic effect with each passage of the roller should be 
considered. 

- The randomness of the surfaces in contact and the loads 
and the material propretles should be incorporated in the 


analytical modelling. 
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APPENDIX A 


Derivation of the Singular Integral 


Because of the singular nature of Kelvin’s solution a cut is 
made in the body to exclude the point P from the region (where 
U. , and T\ . + 00 ), in a Cauchy sense. 

ij iJ 



The line integral going from the surface boundary to the 
singularity point P is cancelled by the integral coming back 
since it can be considered that the same path is being integrated 
over but in opposite direction, therefore Eq . 3.7 will be: 





ds , 





A.l 


where s is the boundary surface of the body and s e is the a 
circle centered at the singularity point P of radius G. 
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Considering the value of the integral at point P as the 
radius of s e goes to zero ( £ 0 ), and noting that : 


ds e = 

e d e 

r ’i = 

n i 

n l = 

- cos 0 

n 2 

- sin 0 


dr/dn= cos 0 + sin0 n^ = ~ 1 

r -i n j • r -j n i ■ ° 
e In e + 0 if e ^ 0 

The values of the integrals around s £ can be evaluated as 
follows : 

Substituting the value of U from Eq . 3.2 


/ t. U. . ds = I 

■J ‘J e J 


6 I 

< (3-4V)]ne 6, . - n.n. [t.edfl A.2 

8ttG(1-v) l 1 J I J 


Placing e outside the Integral: 


{elne|6 t. d0 -e f n.n. t. d © 1 

fl-vl J lj ~ J 1 J J J 


A. 3 


8ttG(1-v)1 


and taking the limit as E -*■ 0, 



U. . ds* - 0 
ij e 


A. A 


By substituting Eq . 3.6 for ^ , the second integral becomes: 
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/ 


u . 
J 


s 


e 




d-2v) 

4tt(1-v) 



+ 


2 

(l-2v) 



ed 6 


A. 5 


(1-2V) 

4tt(1-v) 




_2 

(1-2V) n l n j 


d 0 


« 


Evaluating the function In matrix form: 


(1-2V) 

Air (1-v) 


/* D 

X 1+ T5P 

6 

/ 


cos 2 0 d0 


(1-2V) 


(1-2V) 


sin0 cos0 d0 


Jf 

f'v 


sin0 cos0 dO 


sin 2 0 d0 


A. 6 


Integrating each terra of the matrix, the above equation 


1 


= 4tt(1-v) 


20(l-v)+sin0 cos0 -sin 2 0 

L-sin 2 0 2 0 ( l-V)-sin0cos0 , 


A. 7 


By substituting relations A. 4 and A. 7 in Eq. 
boundary integral equations : 

/ t . U ds = I u . T ds + 

J ij J J iJ 

s s 

r 0 sin0 cos0 sin 2 0 

T 4(l-v) 4 ( 1-v) 


A.l one gets the 


A. 8 


u . 
7T 


sin 2 0 
4 (1-V) 


0 sin0 cos0 

~2 4 (1-v) 



i«y 


If the point P is at the interior of the body, 0 will be equal to 
2tt, if P is a point at a smooth surface, 0 will be equal to TV. 
The integral equation will have the following general form: 


C 


ij 




ds 




ds 


A. 9 


where 


C. . 
ij 


0 sin 0 cosO -sin ? B 

2tT 4tt (1-v) 4tt (1-v) 


♦ 7 r\ 

"Dxn yj 

4 n (1-v) 


n Ql 

2tt 4 it (1-v) I 


A. 10 


where C. . = 6. . for internal points, 
ij ij • 

surface points, with smooth boundary. 


and C. . = 
ij 
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